Preliminary question   If you haven’t yet, read the “Homework set policy” document in the “Prac-ticalities” module on Quercus. Did you read and understand the policy? If you did, write down

“I have read and I understand the homework set policy in this course.”

in your copy. You should answer this question truthfully before continuing with this homework set.

In particular, make sure you re-read the policy about code!

For all questions, if we derived a formula during a lecture or in the lecture notes, you may re-use it without justification other than a simple reference to the notes or textbook.

With this assignment, you should be able to download the following Python scripts:

•  pendulum_solve_ivp.py,

•  find_roots .py, and

•  springpendulum .py.

If any or all are missing, let N.G. know ASAP.

Computational questions can be somewhat open-ended, but this does not mean that you should write more than necessary. Aim for concise answers with clearly labelled and described plots. Nu- merical studies are a bit like experiments; you have to try things and report the most interesting things you observe.  You do not need to hand in every code you use, otherwise it might be a lot. Instead, hand in a few sample codes that represent your work and will help the marker understand why you find what you find.

We highly recommend saving the code you write for each subquestion below with a different name. For example, when you finish question 2.1, save the code as q21 .py, then save a copy of it as q21 .py and alter this version when you move to part 2 of the assignment. Finally, when you upload code, make sure to reference the script in the body of the text. (E.g., “These figures were computed using MyPrettyCode.py.”)

(DO WE?)  We often refer to a certain ddp_routines .py script.  It is part of the material for the lecture “Describing chaos”, which we uploaded on Quercus. The functions contained in this script are more advanced that what we are asking you to plothere, but you can reverse-engineer them to suit your needs. You can do a test run on a few of these functions by running ddp_basic_plots .py, making sure ddp_routines.py is in the same directory first.

If we don’t use the updated value of velocity when calculating position, we are simply using the forward Euler method:

Let’s see what happens when we solve the pendulum using the forward Euler method.

1.1. Write a Python program using the forward Euler method to solve the non-dimensional pen- dulum eq. (1).  Use a time step of 0.05τ and plot the energy and θ (on different plots) as a function of time for 60τ. What do you notice? Does this make sense? Why or why not? Now pick a time step of 0.01τ and do the same, what changes? You can pick your initial conditions (θ(0) and θ(0)) to be whatever you like, but remember that they are in radians and radians/τ . Keep the same initial conditions for each part of this question.

You should see the energy grow withoutbound for any time step that you can take.  While the forward Euler method will produce exact results in the limit Δ t ! 0,for any finite Δ tit is unconditionally unstable for this problem.

1.2.  Repeat part 1.1. with the Euler-Cromer method instead of the forward Euler method. What difference do you see between the results from these two numerical methods?

Figure 1: x (t) from the Rössler system for the first 100 seconds of the systems evolution, solved with the forward-Euler method. Here a = 0.1, b = 0.1, c = 2, x(0) = y(0) = z (0) = 10.0 and Δt = 0.001.

the same figure. In dynamics, systems often have to run for a little while before they collapse onto an attractor. When would you say that this system has approximately collapsed onto its attractor? What is characteristic about this attractor (i.e., describe it)?

2.4.  Change your program to only plot the solution after the system has equilibrated. Now com- pare the results for c = 4 to the results for c = 6.5 and c = 8 (note: if you change c, the equi- libration time may change).  What differences do you notice between these c values?  Pay close attention to how long it takes for the signal to truly repeat itself. Considering z maybe helpful in determining this.

The behaviour you observe as c increases is called period doubling, and it is one pathway by which systems evolve into a chaotic regime.

2.5.  If we continue to increase c we will find that the time it takes for the system to repeat grows longer and longer and eventually reaches infinity (long before c gets to infinity).

We can visualize this by using the local maxima of z as a proxy for the state of the whole system. Modify your program so that it plots all the local maxima of z as a function of c for a range of c ’s between 4 and 15 (see “Steps to modify your program” below). That is, make a plot of z max vs. c (e.g., if there were 6 distinct local maxima for c = 9 you would have 6 points on the z axis above c = 9 on the c axis).  Good results are achieved when you use about 120 values of c between 4 and 15. Run each simulation (one for each c) for about 1300 seconds in order to get enough maxima, and discard the first 250 seconds of each simulation in order to make sure that the solution is on the attractor. This program may take a few minutes to run.

The plot you made is a bifurcation diagram, the regions which are very dense with points correspond to chaotic regions of parameter space. What does this result tell you about the dependence of the system on c?

Steps to modify your program:

•  The function getmaxes given below can be used to find the maxima of z. As usual with Python, you should put this function near the beginning of your program.

def  getmaxes(z,  c):

"""  Finds  the  local  maximas  in  a  given  solution .

Takes  an  array  z  and  a  float  c  and  returns  a  list  of  lists

containing  max list  and  c list .

max list  and  c list  contain  the  values  of  all  the  local  maximas   in  z  and  the  c 's  that  go  with  each  maxima   (the  c 's  will  all  be the  same  for  a  single  getmaxes  call) . """

maxlist  =  []

clist  =  []

a  =  (np.diff(np.sign(np.diff(z)))  <  0) .nonzero()[0]  +  1

for  q  in  range (0 ,  size(a)):

maxlist .append(z[a[q]])

clist .append(c)

return  [maxlist,  clist]

• You will want to call the getmaxes function after you have calculated the values of x, y and z. You can call it using:

#  Call  getmaxes  for  z

tempmaxes  =  getmaxes(z,  c)

maximas .extend(tempmaxes[0])

cs.extend(tempmaxes[1])

•  In order to use this code snippet you will need to define the lists maximas and cs near the beginning of your program with:

#  Define  the  empty  lists

maximas  =  []

cs  =  []

•  Now you will want to create a plot of the maximas vs c.  You can do so by adding the following to the end of your program:

#  Plot  the  c 's  vs  the  maximas

plt.plot(cs, maximas, ' o ' , markersize=1.5)

plt.show()

•  Finally, you will notice that at this point, the program only works for a single c value. You will want to alter your program so you can loop over a range of c values. To do so, you are going to want to define an outer loop to loop over the c values and put your code that updates the x, y, z stuff and calculates the tempmaxes inside the c loop. Just make sure your plotline is not in the c loop (i.e., is not indented).

2.6. Another property of chaotic systems is a sensitive dependence on initial conditions. Go back to the program you wrote for part 2.3. and rewrite your program so that it solves two Rössler systems simultaneously.  You will need two sets of initial conditions, and two sets of time stepping segments. Plot the systems as you did in part 2.3., but plot both x (t)’s on the same subplot, likewise for y(t) and z (t). You can verify things are correct by starting both simula- tions with the same initial conditions, the lines should overlap perfectly.

Add a new subplot which shows the absolute value of the difference between x (t) for the two simulations, if the initial conditions are the same, this difference should be zero.

Note: If you chose the solve_ivp method, you need to force the code to output the solution at fixed times, in order for both arrays of x-values to have the same length and the substraction

to make sense. To do that, create a time array time like we did in class with, for example, the same Δt you used in Q2.2., and add the t_ eval = time option at the end of the solve_ivpcall: solution  =  solve_ivp( ... ,  t_eval=time)

were the ... are a placeholder for the other elements of the call we discussed so far.

Now change one of your x (0) by adding 0.0001 to it but otherwise make all the other param- eters the same as in Q2.3.  (a = 0.1, b = 0.1, y(0) = z (0) = 10, etc.).  Now compare what you see when you set c = 4 and c = 10. What regions of the bifurcation diagram that you found in part 2.5. do these correspond to?

2.7. Write a program to calculate the Lyapunov exponent for the two cases that you looked at in Q2.6. (c = 4 and c = 10) by plotting ln jx1 ¡ x2j for each c, where x1 and x2 are the two slightly different trajectories you plotted in Q2.6..

Remember that in order to calculate the Lyapunov exponent you should only consider the time shortly after your simulation starts. Why is this?

For this question you can eyeball the line of best fitif you like.

3    Springy pendulum [35%]

In this question we will explore the spring pendulum system numerically using Python. This sys- tem consists of a mass on the end of a spring which is free to move in a plane.  No small angle approximation is used in this problem so this system cannot, in general, be solved analytically. For certain parameters this system readily displays chaos, exploring this will be the focus of this question.

This is the other kind of chaos that the lecture notes mention:  Hamiltonian chaos.  There is no forcing, no dissipation, energy is conserved and chaos takes the form of a complete exploration of the accessible phase space rather than collapses onto attractors.

A schematic diagram of this system is shown in figure 2. Here ! is the length the spring stretches to when the weighted system is in equilibrium, !0 is the length the spring is when the unweighted spring is at equilibrium. Finally g and m are the gravity and mass respectively. A coordinate system for this setup can bedefined a number of ways. In this assignment we will use an x-y coordinate system centred at the the equilibrium position of the mass.

Figure 2: A schematic diagram of the spring pendulum system with the mass in the equilibrium position. The coordinate system and important variables are marked.

3.1.  Show that the equations of motion for this system are

Finally, the system may not be constrained by its dynamics to lie on any surface and so the whole phase space will be traversed. In this case the system is in a fully chaotic state.

In class, we could easily identify points on the !ftaxis because that’s the axis along which the time-integration proceeded. Here, on the other hand, the solution does not necessarily have y = 0 exactly at some, or any, time step. What it does have is a set of successive values for y that straddle zero. What you need to do to find a better estimate of the crossing is therefore to interpolate between the values that straddle zero in the right direction.

In the script find_roots .py, we give an example that will find the zeros of x = sin(t), with a time array that does not exactly contain values t = nπ .  To find the zeros of x, it detects the places where two successive values of x change sign in the ascending direction (i.e., it finds the values of i for which xi  < 0 and xi+1 > 0), and interpolates t (meant to approximate arcsin x) between the surrounding values to find an accurate version of the corresponding root of sin(t).

Once you have the times tr for which y(tr ) = 0,y(˙)(tr ) > 0, you can use other interpolations of

x and  at times tr to create the values to display on your Poincaré section.1

Combine find_roots .py and your code from part3.3.  to plot a Poincaré section of (x, ) when y(˙) > 0, y = 0. As during the lecture, you also want to plot them as small dots rather than jointed lines, as in

plt.plot(xPoincare,  dotxPoincare,  ' . ' , markersize=0.5)

You will also want to run them for longer than the previous examples to start seeing some- thing.

Finally, we noticed that the default Runge-Kutta integration scheme implemented by solve_ivp led to inaccurate plots, even for extremely low values of rtol.  While these inaccuracies

were barely perceptible up until now, they add up over long integrations and become vis- ible in Poincaré sections.  Instead, we found that the Radau method gave the best results. (Here for example, we know that in the first case below, Poincaré sections must line up along closed trajectories.) To use the Radau method, simply add it as an option at the end of the solve_ivpcall, like this:

solution  =  solve_ivp( ... ,  rtol=1e-6 , method="Radau")

were “ ...” is a placeholder for the other inputs of the call. You could also consider lowering the relative tolerance, but the Radau method is slow as it is, and we don’t require you to be very precise.

(a)  Set σ = 1/4, E = 1.0 and x0  = −0.75, 0  = 1.0, y0  = 0.20 (you need to figure out y(˙)0  on your own, but hint: checkout part3.5.). Now run your code and examine your Poincaré section.

(b)  Now repeat this but set E = 1.0 and x0 = −0.1, 0 = 1.0, y0 = 0.0.

Can you describe the difference in the dynamics in these 2 cases? 

If we now want to simulate a pendulum that is 1 mm long on a missile that is accelerating at 300 m s¡2, the 0.5 s time step we used previously is nowhere near small enough to resolve the motion whose period is approximately 0.01 s. A time step of 0.5 s is 50 times larger than the period of the pendulum, and so is unable to resolve the motion at all.

The issue is that the time step is dependent on the physics, not on the mathematics.

If we non-dimensionalise the equations then the time step is in units of pendulum period (say 20 time steps per τ) so if the physics of the problem changes, the numerics of the problem don’t.

We are now going to derive a non-dimensional energy for the pendulum. We will start by writing the energy in dimensional form. If we set the origin of the system to beat the equilibrium point of the pendulum my dimensional energy is

In order to non-dimensionalise energy we need to pick a length scale L and a mass scale M. The obvious choices for these are the length and mass of the pendulum. Putting these into eq. (10) (and making sure to non-dimensionalise the E) we get