Each problem counts for 10 points. Full credit for this set = 50 points (do at least all unstarred problems).
1.A stiff system of ODE’s. Use integrate.odeint to compute andplot the solutions u(t), v(t) to the following coupled ODE’s:
over the range 0 < t < 3 with initial conditions u(0) = 2, v(0) = 1. You should see both u(t) and v(t) decay smoothly and slowly towards zero.
Next, solve these equations using second-order Runge-Kutta, plotting u(t), v(t) from RK2 on the same graph as the odeint solutions (with a legend to show what is what, of course). You may want to use something like plt.ylim(-2,4) to fix the y range of the graph.
You should find that you need quite a lot of time points to make theRK2 solution stable, despite the smooth, slow variation of the correct solution. The equations above are an example of a stiff system ofODE’s, which means that there are other solutions to these equations with veryrapid time dependence. For this reason, a very small timestep is required for solution methods like RK2 even though the actual
solution of interest does not have this rapid time dependence. An adaptive step-size method as discussed in Newman Sec. 8.4 will not help with stiff systems; instead more complicated methods are required. The function scipy.integrate.odeint automatically detects stiff systems and switches to an appropriate solver.2
This example describes a pendulum consisting of a mass m connected to a fixed pivot by a rigid rod of length A, and derives the following ODE for the angle θ of the pendulum rod with respect to vertical:
For this problem use m = 1.0 kg, A = 1.0 m, and of course g =m/s2.This problem asks some questions, which you should answer in markdown cells in the Jupyter Notebook you turn in.
(a)Compute and plot θ(t) over the range 0 < t < 5 s with initial conditions at t = 0 of θ = 0, dθ/dt = 1.0 rad/s. Explain your results in terms of the expected period of the pendulum in the linear regime |θ| 1, which you should compute.
(b)Next, change your program to plot θ(t) for 50 different initial conditions all on the same graph (so there will be 50 curves on this plot). For the 50 initial conditions at t = 0, make dθ/dt vary inequal steps from zero to 10 rad/s, while keeping everything else the same as in the previous part. If your variable for the initial value of dθ/dt is called omega0 you can do this by writing
You would put solving the ODE and adding a curve to the figure inside the for loop.
2 I copied this example from Numerical Recipes – The Art of Scientific Computing, Third Edition by Press et al., Sec. 17.5, which has a nice discussion of stiff systems. Systems like this are called “stiff” because they arise in the dynamics systems with stiff constraints between the variables like mechanical linkages. In our example, there is a “stiff linkage” enforcing u(t) 2v(t). Any deviation from this constraint decays away very rapidly, as you can easily test.
What is the exact value of initial dθ/dt that separates the two behaviours? Hint: Energy conservation.
(c)Rather than plot θ or dθ/dt versus t, you can plot θ versus dθ/dt sothese two variables will be the axes of your plot (and there will be no t axis). The plane of dθ/dt versus θ is called phase space, and for every choice of initial conditions the system will trace out a trajectory (curve) in phase space as t varies.3
Plot the same 50 trajectories as computed in the previous part in phase space – you will have 50 curves in the dθ/dt versus θ plane.To see what is going on more clearly, set the θ limits of your plot to ( π, π). You should find that some of the trajectories formclosed curves in phase space, while other trajectories do not. How does this relate to your results from part 2b?
Optional (but not difficult) addition: Fill in all four quadrants of phase space by
i.making your dθ/dt initial values range over ( 10, 10) rad/s (that is, include both negative and positive values),and
ii.addinga second set of curves to your plot going from t = 0 to t = 5 s, that is integrated to negative times. Simply create a second set of time values, e.g. t2 = linspace(0,-5,200) and pass this to odeint.
The code youneed is almost the same as for the “Orbit of a planet” in-class exercise– you just need a different function for the derivatives and you should change the names of the variables to match this problem.
4.Newman8.3: The Lorenz equations (deterministic chaos). To see all the detail in these plots you will want to use a lot of time points. Required extra bit for part (a): Make two graphs of y(t) on the same plot, one starting with the initial conditions Newman lists (0, 1, 0) and one very slightly different initial conditions: (0, 1+10−9, 0).
3 Some picky people will say phase space is really the plane of (mA2)dθ/dt versus θ so that the product of the units of the two axes is distance times momentum which is called “action” (same units as angular momentum and ¯h).
You will see that the two y(t) graphs start out almost identically, as youwould expect. But by the end of the graph they are completely different– an extremely tiny change to the initial conditions has completely changed the results a finite time later.
This is the essence of deterministic chaos. Even though the ODE’s are deterministic (there is nothing random about them), the resultsquickly become unpredictable with any uncertainty in the initial con- ditions no matter how small. The discoverer of these equations, EdwardLorenz, called this the butterfly effect: A butterfly flapping its wings can result in huge changes in the weather some days later. (The sen- sitivity of some deterministic systems like the three-body problem totiny changes in initial conditions had been noted by Henri Poincar´e much earlier, in 1890.)
Plot your x(t), y(t),z(t) data from the previous problem in 3D on x, y, z axes, and make the plot appear in a separate window so you can spin it around with your mouse and see the Lorenz attractor from various angles.4 The code to make this plot is similar to the code for 3D scatter plot in our Plotting and graphics handout, except that instead of scatter(x,y,z) which plots a bunch of points you want ax.plot(x,y,z) which by default will draw a curve connecting the points.
6.*Hit a target with the cannonball. This continues problem 3 above.Use a bisection search to adjust the speed at which the can- nonball is fired to precisely hit a target some specified distance See Newman Sec. 8.6 Boundary value problems to get started (Newman calls a bisection search a “binary search”). Note: This is rather more challenging than most problems we have done. I was able to make it work – if you try this, you may want to get some hints from me.
4For the Jupyter notebook you turn in, you should put a 3D plot in the notebook –even better, a couple of 3D plots from different angles.
