6810: 1094 Session 8
Online handout: plots of damped oscillations;
filename_test.cpp, diffeq_pendulum.cpp, GnuplotPipe class
Strings and Things
The filename_test.cpp code has examples of the use and manipulation of
C++ strings, including building filenames the way we do stream output.
Be careful NOT to put << endl when creating filenames.
- Using make_filename_test, compile and link filename_test.cpp and
run it. Look at the output files and the printout of the
code to see how it works.
- Modify the code so that there is a loop running from 0 to 3 with
index variable j. For each j, open a file with a name that includes
the current value of j.
Write "This is file j", where "j" here is the
current value, into each file and then close it. Did you succeed?
- Modify the code to input a double named alpha
and open a filename with 3
digits of alpha as part of the name. (E.g., something like
pendulum_alpha5.22_plot.dat if alpha = 5.21934.) Output something
appropriate to the file. Did it work?
Upgrades from the diffeq_oscillation to diffeq_pendulum code
- There are three new menu items: plot_start, plot_end, and
equation is still solved from t_start to t_end, but results are only
printed out from plot_start to plot_end. Initially these are the same time
intervals, but you can use plot_start to exclude a transient
So if the system settles down to periodic behavior at t=20, setting
plot_start=20 means that 0 < t < 20 is not plotted, which makes the
phase-space plots much easier to interpret.
- We've also incorporated code to do real-time plotting in gnuplot directly from
C++ programs. We have made a class to do this but it is
and documentation needs work, and it probably has bugs!
Look at the GnuplotPipe.h printout and the GnuplotPipe.cpp file
to get an idea how it works.
Gnuplot_delay sets the time in milliseconds between plotted points.
Damped (Undriven) Pendulum
The pendulum modeled here has the analog of the
viscous damping: Ff = −b*v, where v(t) is the velocity, that
was used in session 7. The damping parameter is called alpha here.
- Use make_diffeq_pendulum to compile and link diffeq_pendulum.cpp.
Run it while taking a look at the printout of the
It should look a lot like
diffeq_oscillations.cpp, with different parameter names. Run it with
the default parameters, noting the real-time phase-space plot. There
is also an output file diffeq_pendulum.dat.
- Modify the code so that the output file includes two digits of
the variable alpha in the name. Did you succeed?
- Generate the analogs of the four phase-space plots on the handout
but with pendulum variables and initial conditions theta_dot0=0 (at
rest) and theta0 such that you are in the simple harmonic oscillator
regime (note that theta is in radians). Set f_ext=0 (no external
driving force) and then do four runs with four values of alpha
corresponding to undamped, underdamped, critically damped, and
overdamped (convert from the conditions on b discussed in the
background notes). What values of theta0 and alpha did
Damped, Driven Pendulum
This is a quick exercise to look at transients.
- Restart the program so that we use the defaults. There is both
damping and an external driving force, with frequency w_ext = 0.689.
The initial plot is from t=0 to t=100. Run it.
The green points are
plotted once every period of the external force. What
good are they?
- Note that it seems to settle down to a periodic orbit after a
while. Plot ("by hand" with gnuplot) theta vs. t from the output file
diffeq_pendulum.dat and see how long it takes to become periodic.
- Run the code again with "plot_start" set to the time you just
Have you gotten rid of the transients? What is the frequency
of the asymptotic theta(t)?
Looking for Chaos
Now we want to explore more of the parameter space and look at different
structures. In Section f of the Session 7 notes
there is a list of characteristic structures
that can be found in phase space, with sample pictures in Figure 1.
- In phase space, a fixed point is a (zero-dimensional) point that
"attracts" the time-development of a system. By this we mean that
many (or all) initial conditions end up at the same point in phase
space. The clearest case is a damped, undriven system like a pendulum,
which ends up at theta=0 and zero angular velocity no matter how it
starts. If the steady-state trajectory in phase space is a closed
(one-dimensional) curve, then we call it a limit cycle.
- Try some prescribed values for the pendulum.
You will need to adjust "plot_start" and extend the plot time
(increase "t_end" and "plot_end").
Try the first three combinations in this table:
|period-1 limit cycle||0.0||0.0||0.689||0.8||0.0|
Can you tell how many "periods" the limit cycles have from the
graphs? How might you identify whether a function of time f(t) is
built from one, two, three, ... frequencies?
One characteristic of chaos is an "exponential sensitivity to initial
conditions." For the last combination, vary the initial conditions
(e.g., change x0 by 0.01 or 0.001); what happens?
Armadillo linear algebra library (Do this on Linux!)
Here we'll use the Armadillo library as an example of how to install a personal copy
of a library and set up a makefile appropriately to use it. We'll try some basic examples
and see how to use it as an alternative to GSL in the Hamiltonian class.
- Installing Armadillo. Follow these instructions but ask if something goes wrong.
- Create a subdirectory called my_armadillo in your 6810 directory on Linux.
- Download from
the latest Armadillo tarball (named something like armadillo-#.###.#.tar.gz)
into my_armadillo and unpack it with tar xfvz armadillo-#.###.#.tar.gz
(substitute for the #'s appropriately).
- Go into the resulting directory (named armadillo-#.###.#) and look at
README.txt. Give the command pwd to find the full path to this directory;
we'll need to use it several times so I'll call it <path-to-Armadillo>.
When I tested it, my <path-to-Armadillo> was:
- Now type (with returns after each command!) cmake . (note the period), then make, then
make install DESTDIR=. (note the period again). This first creates a makefile with
cmake, then compiles the library, than installs it in the same directory (we could put it
somewhere else by replacing the last period by a path to another directory).
Did this seem to work?
- Testing Armadillo. Go to the session_08 directory and
edit section 3. of make_armadillo_tests to prepare it for your installation.
You can browse the online documentation
for other examples of using the Armadillo library.
- We need to link to three libraries. After LIBS= add
-larmadillo -llapack -lblas
- We need to set the path to the Armadillo library. After LDFLAGS= add
- We need to set the path to the Armadillo headers. After CFLAGS= -g -O2
- Now try to compile and link armadillo_tests. It should work to create the .x file,
but fail when you
try to run armadillo_tests.x. The problem is that you need to set the path to the Armadillo
library. At the terminal prompt, type:
setenv LD_LIBRARY_PATH <path-to-Armadillo>
if you are using tcsh. If you are using bash (type echo $SHELL to check), type:
Now it should work!
- Look at the code and play a bit (e.g., modify the examples). Then add code to solve:
x1 + 2x2 + 3x3 = 2
2x1 − 3x2 + 4x3 = 5
x1 + 4x2 − 2x3 = −2
as converted to the matrix problem
Ax = b.
What is the result for the solution vector x?
- Armadillo for the Hamiltonian class. An alternative version of eigen_tridiagonal_class.cpp,
called eigen_tridiagonal_class_armadillo.cpp, uses Armadillo as an alternative to GSL.
- The diff file1 file2 command can be used to show all of the line-by-line
differences between two files file1 and file2.
for more details on diff. Type man diff to see all of the (many) options.]
Use diff to compare
the two eigen_tridiagonal cpp files (both are in session08.zip). Why are there
so few differences? What has been hidden?
- Look at the ArmadilloHamilonian header and cpp files, looking up the functions
in the Armadillo documentation (see the comments). Which elements of the class do you
- Edit make_eigen_tridiagonal_class_armadillo to use your personal Armadillo library
and compile, link, and run eigen_tridiagonal_class_armadillo. Compile, link, and run the
GSL version to create another .dat file and use eigen_tridiagonal_comparison.plt to
compare the outputs. Are they the same? Is this great, or what?
6810: 1094 Session 8.
Last modified: 08:49 am, February 20, 2017.