MECH 202 Jupyter Notebook Written Report – Van der Pol
fouMECH 202
Jupyter Notebook Written Report
Van der Pol Oscillator
Due date: Sunday of Week 12
Grading & Weight: This assignment is out of 100 marks, as further specified in the rubric at the
end of this document, and is worth 8.5% of your overall final grade in the course. There are 12
bonus marks available.
Late Penalty: Late submissions will be penalized at 10% each day for up to 5 days in which
case a grade of zero will be given.
1. Overview
For this assignment, you will solve a non-linear second-order ordinary differential equation using
numerical methods. As explained in class, this project will be completed over two weeks:
Part 1 (Week 11): - Simple oscillator
- Damped Oscillator
- rewrite Van der Pol as system of 1st order ODEs
- Solve for three particular cases
Part 2 (Week 12)
-Plot velocity-position diagrams
-Driven Van der Pol oscillator
-Bifurcation diagrams
You will submit a report in the form of a Jupyter notebook.
This assignment directly aligns with the following Course Learning Outcomes (CLO):
CLO (5): Solve ordinary differential equations numerically, coding in Python.
1.1.Time for completion
This assignment will take approximately twelve (12) hours to complete. You will begin working
on this assignment in the two-hour tutorials of week 11 and 12, and finish on your own time
before the deadline.
MECH 202 Jupyter Notebook Written Report – Van der Pol
2. Instructions
For this assignment, you will implement the numerical integration techniques that we saw in the
online content and in-class. You will apply this on various 2nd order ODEs. You will submit a
report in the form of a Jupyter notebook. Your Jupyter notebook must contain the following
information:
Introduction
o In your introduction, include pertinent theoretical background. This can be
completed in about three paragraph (200-300 words). Include references as you
judge pertinent. You might want to refer to some of those presented in class.
Completion of the tasks and discussion
o Include code required to correctly complete the task
o Include comments to explain your thought process and code logic
o Discuss the results of what was found and answer the questions in the task list
Conclusion
o Discuss what you learned, what surprised you, and how you could extend this
work. This can be completed in about 250 words.
Optional reference list
o Include a references if sources are used
o This can be included at the end of the intro, or at the end of the document. This
esthetic choice is up to you!
Format
o See the report template in onQ to format your report.
Part 1 (WEEK 11)
Part 1: Simple Oscillator:
The simple oscillator can be written as:
𝑞 + 𝑘𝑞 = 0
In mechanical engineering, it is often seen as:
𝑥 + 𝜔02𝑥 = 0
Steps:
(for this part, use your own ODE solver; do not use the ones already implemented in Python unless stated
otherwise)
1. Rewrite this second-order ODE as two first-order ODEs (using Jupyter Markdown) (3 marks)
MECH 202 Jupyter Notebook Written Report – Van der Pol
2. Set the frequency at π/8 (rad/s). Set the initial position at x=1 (arbitrary units), and velocity
at v = 0. Solve the system for its first 180 seconds. Use a dt of 0.1s. (Present a graph with three
curves, each corresponding to an ODE scheme). Use:
a. Euler’s method (3 marks)
b. Explicit Runge-Kutta 4th order (3 marks)
c. SciPy’s solver – make sure to use a version that handles stiff equations (3 marks)
Part 2: Simple Damped Oscillator
The simple damped oscillator can be written as:
𝑞 + 2𝜁𝑞 + 𝑞 = 0
It is often seen as:
𝑥 + 2𝜁𝜔0𝑥 + 𝜔02𝑥 = 0
with
𝜁 = 𝑐 2√𝑚𝑘
and
𝜔0 = √𝑘𝑚
Steps:
(for this part, use your own ODE solver methods; do not use the ones already implemented in Python
unless stated otherwise)
1. Rewrite this second-order ODE as two first-order ODEs (using Jupyter Markdown) (3 marks)
2. Set the frequency at π /8 (rad/s). Set the damping ratio at 1/16. Set the initial position at
x=1 (arbitrary units), and velocity at v = 0. Solve the system for its first 180 seconds. Use a dt of
0.1. (Present a graph with three curves, each corresponding to an ODE scheme). Use:
a. Euler’s method (3 marks)
b. Explicit Runge-Kutta 4th order (3 marks)
c. SciPy’s solver – make sure to use a version that handles stiff equations (3 marks)
3. Set the frequency at π/8 (rad/s). Set the damping ratio at 4. Set the initial position at x=1
(arbitrary units), and velocity at v = 0. Solve the system for its first 180 seconds. For each solver,
try dt of 0.1. (Present one graph with three curves, each corresponding to an ODE scheme). Use:
a. Euler’s method (3 marks)
b. Explicit Runge-Kutta 4th order (3 marks)
c. SciPy’s solver – make sure to use a version that handles stiff equations (3 marks)
d. What do you observe? What can you say about the stability of Euler’s and RK-45? You’ll
notice that stability issues arise; in other words, this is a stiff mechanical system. In
general, the term “stiff equation” is used to describe these types of rapidly-varying
systems (including equations linked to fluids, neurobiology and electronics). (3 marks)
4. Set the frequency at π/8 (rad/s). Set the damping ratio at -1/16 (i.e. negative damping).
Set the initial position at x=1 (arbitrary units), and velocity at v = 0. Solve the system for its first
180 seconds. Use a dt of 0.1 (Plot all of the results on the same graph). Use:
a. Euler’s method (3 marks)
MECH 202 Jupyter Notebook Written Report – Van der Pol
b. Explicit Runge-Kutta 4th order (3 marks)
c. SciPy’s solver – make sure to use a version that handles stiff equations (3 marks)
Part 3: Van der Pol oscillator
Now we’ll implement the Relaxation-Oscillation equation from the paper (Van der Pol, B., 1926.
LXXXVIII. On “relaxation-oscillations”. The London, Edinburgh, and Dublin Philosophical
Magazine and Journal of Science, 2(11), pp.978-992.) which can be found at
https://www.tandfonline.com/doi/abs/10.1080/14786442608564127?journalCode=tphm18.
The idea is to have a damping which is negative when the oscillator is near its equilibrium
position, but positive when it is further away.
The relaxation-oscillation equation is given by:
𝑥 − 𝜁(1 − 𝑥2)𝑥 + 𝑥 = 0
Steps:
(for this part, use your own ODE solver methods; do not use the ones already implemented in Python
unless stated otherwise)
1. Rewrite this second-order ODE as two first-order ODEs (using Jupyter Markdown) (2 marks)
2. Set the initial position at x=0 (arbitrary units), and velocity at v = 0.02. Solve the system for its first
180 seconds. . Use a dt of [0.01] and [0.5]. Use Euler’s method. For each , plot one graph of
position as a function of time. :
a. Set at 0.1 (3 marks).
b. Set at 1.0 (3 marks).
c. Set at 10 (3 marks).
You’ll likely notice that using a dt of [0.5] doesn’t work in all situations. Why is this?
3. Set the initial position at x=0 (arbitrary units), and velocity at v = 0.02. Solve the system for its first
180 seconds. . Use a dt of [0.01] and [0.5]. Use explicit Runge-Kutta 4th order. For each , plot
one graph of position as a function of time :
a. Set at 0.1 (3 marks).
b. Set at 1.0 (3 marks).
c. Set at 10 (3 marks).
You’ll notice that using a dt of [0.5] doesn’t work in all situations. Why is this?
4. Set the initial position at x=0 (arbitrary units), and velocity at v = 0.02. Solve the system for its first
180 seconds. . Use a dt of [0.5]. Use SciPy’s solver – make sure to use a version that handles stiff
equations. For each , plot one graph of position as a function of time. :
MECH 202 Jupyter Notebook Written Report – Van der Pol
a. Set at 0.1 (3 marks).
b. Set at 1.0 (3 marks).
c. Set at 10 (3 marks).
Part 2 (Week 12)
Part 4: Van der Pol oscillator: position-velocity diagram
In this part, we’ll look at the cyclic behavior of the systems by drawing diagrams of velocity as a
function of position.
Steps
1- Set the initial position at x=0 (arbitrary units), and velocity at v = 0.02. Solve the system for its first
180 seconds. . Use a dt of [0.1]. Use SciPy’s solver – make sure to use a version that handles stiff
equations. For each , plot one graph of velocity as a function of position :
a. Set at 0.1 (3 marks).
b. Set at 1.0 (3 marks).
c. Set at 10 (2 marks).
d. In all three cases, the system starts at (x = 0, v = 0.02). How would you describe its
behavior thereafter? Does it eventually reach a steady-state? (2 marks)
Part 5: Forced Van der Pol oscillator (BONUS)
In the previous scenarios we examined the natural oscillatory behaviour of the Van der Pol Oscillator
at various damping ratios (ζ) with the initial state of (x = 0.0, 𝑣=0.02).
𝑥 − 𝜁(1 − 𝑥2)𝑥 + 𝑥 = 0
In the following section we can examine the scenario where we apply an external force to the
oscillator. In this scenario we expose the system to a forcing function.
𝑥 − 𝜁(1 − 𝑥2)𝑥 + 𝑥 = 𝐹
where F is typically of the form 𝐴 ∗ 𝑐𝑜𝑠(𝜔𝑡 + 𝜙)
Rearranging the equation, we get:
𝑥 − 𝜁(1 − 𝑥2)𝑥 + 𝑥 − 𝐴 ∗ 𝑐𝑜𝑠(𝜔𝑡 + 𝜙) = 0
Steps:
(you can use a built-in Python solver)
1- Rewrite this second-order ODE as two first-order ODEs (using Jupyter Markdown) (2 marks)
MECH 202 Jupyter Notebook Written Report – Van der Pol
2- Set the initial position at x=0 (arbitrary units), and velocity at v = 0. Set at 8.53. Set A at 1.2. Set
at 1/5 Pi. Solve the system for its first 180 seconds. . Use a dt of 0.1. Use a SciPy solver.
a. Plot position as a function of time. (2 marks)
b. Plot velocity as a function of velocity. (2 marks)
c. Does the velocity-position diagram follow the same behavior than in the unforced
scenario? Does it reach a steady orbit? (2 marks)
Part 6: Bifurcation diagram (BONUS)
In the previous sections we saw differing behaviour of the Van der Pol oscillator depending on the values
of the parameters over time, but it may be of interest to investigate the behaviour of the Van der Pol
oscillator as a function of a parameter. In this section we will examine the behaviour of the Van der Pol
oscillator as a function of the forcing frequency (𝜔). Since we’ll be examining the behaviour of the
oscillator as a function of 𝜔, we’ll need to adjust the way we view our output variable (position, x(t)) so
that we can create the plot with respect to 𝜔.
We know that the oscillator is periodic and “repeats” itself at steady state, we will view the steady
state behaviour of the oscillator and its response to the changing parameter values. These types
of plots are known as a Pointcaré map.
Steps:
(you can use a built-in Python solver) Pointcaré map of 𝜔
1. Set the initial position at x = 0 (arbitrary units), and velocity at v =0. Set at 3 and A at 5. Vary
(rad/s) from 0.1 to 1 with increments of 0.1 and from 1 to 7 with increments of 0.01. Define
the sampling period as 𝑑𝑡 = 2𝜋⁄𝜔 (note the unit conversion). In other words, store the value of
the position every dt. The idea is that—if the oscillator has linear response—it wil oscillate at a
frequency that is a multiple of the forcing frequency. Otherwise, it will react non-linearly and even
exhibit chaotic behaviour. Set the maximum observation time as 𝑡𝑚𝑎𝑥 = 1000 ∗ 𝑑𝑡 and the
minimum observation time as 𝑡𝑚𝑖𝑛 = 𝑡𝑚𝑎𝑥 − 50 ∗ 𝑑𝑡 such that you obtain 50 samples for each
forcing frequency. Using a SciPy solver:
a. (Bifurcation plot) - Plot the set of positions as a function of . (2 marks)
b. Comment on observations from your Bifurcation plot. Can you identify when the oscillator
exhibits chaotic behaviour? (2 marks)
MECH 202 Jupyter Notebook Written Report – Van der Pol
Your assignment will be evaluated using the following criterion:
Criteria Mastery (A+) High Quality (A) Developing (B) Marginal (C) Not