Lecture 7

Continuous System Simulation

M&S
Continuous systems simulation
▶ Time is treated as a continuous variable that drives the
simulation
▶ Model is based upon differential equations, which describe how
systems evolves over time, and how it responds to changes in
input variables
▶ In this course, we will mostly deal with Ordinary Differential
Equations (ODE), though Partial Differential Equations are
used in some cases

2 / 77
Continuous system simulation: Variables
▶ Three set of variables:
▶ State variables;
▶ Input variables; and
▶ Output variables.
▶ There is no overlap between input and state variables
▶ Input variable are not controlled by the simulation
▶ The angle of the steering wheel in the car simulator example
▶ Output variable are things that we observe
▶ The speed of the vehicle in the car simulator example
▶ State Variables are controlled by the differential equations
▶ The speed of the vehicle in the car simulator example
▶ The location of the vehicle in the car simulator example

3 / 77
Ordinary Differential Equations (ODEs)
An ODE consists of the following ingredients:
▶ An independent variable (usually “time” t that derivatives are
taken with respect to
▶ A dependent variable, i.e. function of the independent variable,
e.g. x = x(t) “the variable x which is a function of t”.
▶ A multi-variable function F that describes a relationship
between the derivatives of the dependent variable (taken with
respect to the independent variable)

Putting it all together, we get

!
dx d2 x d(n−1) x dn x
F t, x, , 2 , · · · , (n−1) =
dt dt dt dtn

4 / 77
ODE: comments
dtn denotes derivative of x w.r.t. t
▶ dx
▶ ddtnx denotes n-th derivative of x w.r.t. t
▶ x′ = dxdt
▶ Dx = dx dt
▶ Order of a differential equation is the highest order of
derivative in that equation

Examples
▶ mx′′ = F (order is 2)
▶ x′ + 32x′′ + x′′′ = 0 (order is 3)
▶ x′ + 34x = 32 (order is 1)

5 / 77
Solving differential equations
▶ Solution is the dependent variable x = x(t) that satisfies the
equation
▶ Key idea: integration

Example: x′′ = 2
1. Integrate once: x′ = 2t + C1
2. Integrate again: x = t2 + tC1 + C2 (Solution)

How to solve for C1 and C2 , also called, constants of integration?

▶ Use initial or boundary conditions.
▶ Using initial conditions x′ (0) = 3 and x(0) = 2, we get C1 =3
and C2 = 2. The solution is x(t) = t2 + 3t + 2.

6 / 77
Aside: Constant of Integration
Notice that
df (x)
= 3x2
dx
for both when f (x) = x3 or f (x)x3 + C.
This suggests constant C disappears in the process of differentiation.
Therefore, when we integrate we add the constant C for the sake of
completeness.
!
x3
Z
2
3x dx = 3 + C = x3 + C
3

Furthermore, we will need other information to find the true value

of C. Note also that there is nothing preventing C = 0.

7 / 77
Nth order ODEs
▶ General solution to an nth order ODE will contain n constants
of integration
▶ We need n more equations
▶ Use initial or boundary conditions to get these equations to
solve for the constants of integration
▶ Initial conditions
▶ The values of x(t) and its first n − 1 derivatives for a particular
value of t
▶ If such values are only available at the end, run time backwards
to convert problem to initial conditions
▶ Boundary conditions
▶ The values of x(t) and its derivatives for two different values of
t

8 / 77
Nth order ODE reducibility
Any explicit differential equation of order n
 
x(n) = F t, x, x′ , · · · , x(n−1)

can be written as a system of n first-order differential equations by

defining a new family of unknown functions

x(i−1) = xi

Notice the abuse of notation. Here x(k) denote the k-th derivative
of x w.r.t. t.

9 / 77
Nth order ODE reducibility
We can then represent the following n-th order ODE
 
x(n) = F t, x, x′ , · · · , x(n−1)

with n first-order ODEs as follows

x′1 = x2
x′2 = x3
..
.
x′(n−1) = xn
x′n = F (t, x1 , x2 , · · · , xn )

10 / 77
Nth order ODE reducibility
Example
The following 2nd order ODE

d2 x
= −g
dt2
can be reduced to
dx
=v
dt
dv
= −g
dt
We introduced a new variable v.

11 / 77
First order ODEs
▶ All of our simulations only involve first order ODEs
▶ What about models that involve higher order ODEs?
▶ E.g., the equation of motion for particle is modelled by a second
order differential equation
▶ Applies reducibility, i.e., replace a higher order differential
equation by a system of first order differential equations
▶ We can replace an nth order ODE with n first order ODE
▶ Advantages of first order ODEs
▶ First order equations are much easier to solve numerically
▶ Very few numerical solvers available for higher order equations

12 / 77
Equation of motion
Newton’s second law of motion:
“The acceleration of an object as produced by a net force
is directly proportional to the magnitude of the net force,
in the same direction as the net force, and inversely pro-
portional to the mass of the object.”

Mathematically:
a ∝ F and a ∝ m 1
, and combining the two we get F = ma. Recall
2
a = dt2 , so F = ma is a second order equation.
d x

F = ma
d2 x
=⇒ F = m
dt2

13 / 77
Equation of motion
We introduce a new variable velocity v = dt ,
dx
and get the following
first order system of equations

dx
v=
dt
dv
F =m
dt

14 / 77
Solving ODEs numerically
General idea: given a solution x(t) at time t = t0 , incrementally
step forward in time to find x(t + ∆t)
Example: lets consider the equation of motion

dv F
 
F =m =⇒ ∆v = (∆t)
dt m
dx
v= =⇒ ∆x = (v)(∆t)
dt
We can use ∆v and ∆x to update the current values of v and x,
respectively.

15 / 77
Solving ODEs numerically: an example
What is the value of x at time t = 3?
Model

v(t + ∆t) = v(t) + ∆v = v(t) + (F/m)(∆t)

x(t + ∆t) = x(t) + ∆x = x(t) + (v)(∆t)
Setup
▶ m = 1, F = 1 (other quantities used in the simulation)
▶ v(0) = 0, x(0) = 0 (initial state), ∆t = 1 (time step)
Solution
▶ v(1) = 1
▶ x(1) = 1
▶ v(2) = 2
▶ x(2) = 3
▶ v(3) = 3
▶ x(3) = 6
16 / 77
Solving ODEs numerically: practical considerations
▶ Have to choose ∆t carefully
▶ ∆t too small; the simulation can become very slow
▶ ∆t too large; the simulation can become very inaccurate
▶ Advanced techniques can change ∆t when solving equations to
maintain acceptable accuracy and speed

▶ ∆t determines the exact points in time for which we have the

solution
▶ What if we want solution at other points in time?
▶ This places constraints on how we solve the equations

▶ Often times ∆t used for solving ODEs is much smaller than the
one used to update the display

17 / 77
Choice of ∆T
▶ Molecular activity (fraction of a millisecond)
▶ Evolution of an ecosystem (months or years)
▶ Galaxy formation (millions or billions of years)

18 / 77
Simulation loop
▶ Advance the simulation
▶ Display current results
▶ Get the user response
▶ Repeat

19 / 77
Mass spring system
Hook’s law (1676) states, “the
extension is proportional to the
force.”

Mathematically, F = −kx, where k

is the spring constant and x is the
displacement of the spring from rest
position under the application of
force F .

21 / 77
Mass spring system
Step 1: construct a model that will describe the motion of
the mass over time
Hook’s law: F = −kx
Newton’s Second Law of Motion: F = ma
Combining the two we get

ma = −kx
dx2
=⇒ m = −kx
dt2

22 / 77
Mass spring system
Step2: find a way to solve the model numerically
2
Convert the second order m dx
dt2
= −kx to a system of first order
equations
dx
=v
dt
dv
m = −kx
dt
And make the update rules

x(t + ∆t) = x(t) + v(t)∆t

k
v(t + ∆t) = v(t) − x(t)∆t
m

23 / 77
Mass spring system
▶ Set values for mass m and spring constant k
▶ Set the initial conditions
▶ Values for x and v at start time t0
▶ Run the simulation loop
1. Update t to t + ∆t
2. Update values for x and v using the update ruls
3. Display results or save them to file for plotting
4. Repeat steps 1 to 4

We just simulated a Simple Harmonic Oscillator

eeeeeeeeee
eiiiiiiiiiiiiiiieeeeeeeeeeeeeeeeeeeeeeeeeeeeee
wwwwwww
sdsdsdsd
wwwwwwwwwwwwwwww
sdssdsdd
wwwwwww
Code: 1d-mass-spring

24 / 77
Mass spring system
# Mass-Spring system
class Mass:
def __init__(self):
self.x = 5
self.vx = 0
self.k = 1
self.dt = 0.1
self.t = 0
self.m = 1.0

def update(self):
self.x += (self.vx * self.dt)
self.vx += (- self.k * self.x * self.dt / self.m)
self.t += self.dt

25 / 77
Mass Spring Damper
The mass experiences a damping force that is proportional to its
current velocity

Mathematically
F = −kx − cv
where c is the damping constant

Code: modify 1d-mass-spring to add damping effect

26 / 77
Bouncing ball
+!
Assumption 1: We simplify
the problem by treating the &
ball as a particle
From Newton’s Second Law of
Motion ℎ

dx2
F =m !=0
dt2
where x is the height of the
ball from the ground and m is
the mass of the ball.

27 / 77
Bouncing ball
Assumption 2 Gravity is the +"
only force acting upon this
'
ball then !

F = −mg
ℎ
Putting it together we get

dx2
= −g "=0
dt2

28 / 77
Bouncing ball
Data collection +"
▶ We need to know the
' !
value of g. For our
purposes, we use
g = 9.8m/s2
ℎ
▶ By using different g we
can simulate bouncing
ball on different planets
"=0

29 / 77
Bouncing ball

Did you notice something peculiar with this plot?

Code: ball-floor turn off RK4

30 / 77
Bouncing ball
▶ The ball goes higher with each bounce, which is unexpected.
▶ The error doesn’t go away even if we make timestep really
small. It does, however, minimizes the effect.
▶ It seems we are imparting energy to the ball with each bounce.
This breaks the the law of conservation of energy, which states
that “the total energy of an isolated system remains constant.”

31 / 77
Bouncing ball total energy
Total energy of the ball is the sum of its kinetic and potential
energies.
Kinetic energy = 21 mv 2
Potential energy = mgy

This behavior is due to incorrect assumptions of Euler method. 32 / 77

Bouncing Ball
Total energy is conserved when using Runga-Kutta or RK4 solver.

Code: ball-floor turn on RK4

33 / 77
Euler method
▶ A numerical solver for first order ODEs
▶ First order numerical procedure for solving ODEs (initial value
problems)
▶ It is an explicit method
▶ Calculates the state of the system at a later time given its
current state by using the update equations
▶ y(t + ∆t) = F (y(t))

Aside: implicit methods

▶ Calculates the state of the system at a later time given by
solving an equation that includes both the future state and the
current state
▶ G(y(t + ∆t), y(t)) = 0

34 / 77
Euler method
▶ Numerically unstable
▶ Adversally effects accuracy
▶ Exhibits error growth over time
▶ Error is proportional to ∆t
▶ Particularly unsuited for stiff equations
▶ Equations containing terms that lead to rapid changes
▶ E.g., a mass spring system with large spring constant
▶ Use extremely small time steps
▶ Infeasible in practice

35 / 77
Runga-Kutta method
▶ An other numerical solver for first order ODEs
▶ An alternate to Euler method
▶ A family of explicit and implicit methods
▶ Often RK4 is used
▶ Error is proportional to ∆t4
▶ Makes a huge difference for small values of ∆t

Takeaway: whenever possible use RK4 method

36 / 77
Numerical solvers in Python
from scipy.integrate import ode
def f(self, t, y, arg1):
"""Solves y' = f(t, y)
Arguments:
- y is the state of the system. In our case
y[0] is the position and y[1] is the velocity.
- arg1 is 9.8, as set by set_f_params() method.

Returns vector dy/dt. In our case, dx/dt = v and

dv/dt = -g.
"""
return [y[1], -arg1]

r = ode(f).set_integrator('dop853')
r.set_initial_value([y0, vy0], t0)
r.set_f_params(9.8)
r.integrate(dt)
print r.t, r.y
37 / 77
Bouncing ball: takeaways
▶ Exploit your knowledge of physics to determine if simulation is
behaving as expected
▶ Use several strategies
▶ Compare outputs of several strategies
▶ If outputs differ, you must have a way to explain the differences
▶ If outputs are the same, the simulation may be correct

38 / 77
Discussion
Q. Why does Euler method performing so poorly for our bouncing
ball example?
A. Euler method assumes that the acceleration remains constant
between two time steps. Notice that this assumption is generally
false, but especially so when the ball “hits” the ground at x = 0.
The velocity is flipped, changing the sign of the derivative and
causing a discontinuity.
RK4 method is much better at handling discontinuities (as long as
there aren’t too many of these).
This is why RK4 is able to get good results even for large time steps.

39 / 77
Bouncing ball

For this simulation, the floor sits at height 0. The ball pierces
through the floor, which is incorrect.

Code: ball-floor increase timestep to see the ball penetrating the floor

40 / 77
Bouncing ball
Need a better way to detect collisions with the floor
Scheme 1
▶ Use smaller time steps
▶ The ball will travel less distance between two time steps, and
there is a greater chance of catching the collision instant
▶ In any case, the ball will penetrate less into the floor

Scheme 2
▶ Try to find the exact time of collision using x = vt relationship
▶ Adjust time step accordingly

41 / 77
Bouncing ball
Collision detection
1. Approximate time to collide tc = x
v
2. Set x = 0 and v = −v(t + tc )
▶ Flip v to indicate that the ball is now going back up again

Problem
v is larger than had we calculated tc exactly right (that’s because
the particle is under constant acceleration). Consequently energy is
not conserved.

42 / 77
Bouncing ball
▶ Use the Law of Conservation of Energy to compute the velocity
of the ball when it touches ground.
▶ The ball was released at height h. We know the total energy of
the system, which is mgh. At the start the kinetic energy is 0.
▶ When the ball touches the ground, its potential energy reduces
to 0. Since the total energy remains the same, all of its energy
is now kinetic energy.
1
mv 2 = mgh
2 p
v = 2gh
√
i.e., set x = 0 and v = − 2gh at collision time.

43 / 77