Physics 281: Computational Physics
Homework Set 7
Due in Moodle Thurs 4/11/191
Note: These problems are all about solving ODE’s. The rst problem asks you to compare the performance of scipy.integrate.odeint and second-order Runge-Kutta.
For all the remaining problems, you can just use scipy.integrate.odeint to numerically solve the ODE’s.
Homework problems
Each problem counts for 10 points. Full credit for this set = 50 points (do at least all unstarred problems).
1.A sti system of ODE’s. Use scipy.integrate.odeint to compute and plot the solutions u(t), v(t) to the following coupled ODE’s:
du = 998u + 1998v; (1)
dt
dv = 999u 1999v (2)
dt
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 x the y range of the graph.
You should nd that you need quite a lot of time points to make the RK2 solution stable, despite the smooth, slow variation of the correct solution. The equations above are an example of a sti system of ODE’s, which means that there are other solutions to these equations with very rapid time dependence. For this reason, a very small time step is required for solution methods like RK2 even though the actual
12019-04-02dr P281 HW7.tex c 2019 Donald Candela
1
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 sti systems; instead more complicated methods are required. The function scipy.integrate.odeint automatically detects sti systems and switches to an appropriate solver.2
2.The nonlinear pendulum. Read Example 8.6 on pp 349-350 of Newman. This example describes a pendulum consisting of a mass m connected to a xed pivot by a rigid rod of length , and derives the following ODE for the angle of the pendulum rod with respect to vertical: d2 g dt2 = sin :
For this problem use m = 1:0 kg, ` = 1:0 m, and of course g =
9:8 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 j j 1, which you should compute.
(b)Next, change your program to plot (t) for 50 di erent 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 in equal 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
for omega0 in np.linspace(0,10,50):
You would put solving the ODE and adding a curve to the gure inside the for loop.
2I copied this example from Numerical Recipes – The Art of Scienti c Computing, Third Edition by Press et al., Sec. 17.5, which has a nice discussion of sti systems. Systems like this are called sti ” because they arise in the dynamics systems with sti constraints between the variables like mechanical linkages. In our example, there is a sti
linkage” enforcing u(t) 2v(t). Any deviation from this constraint decays away very rapidly, as you can easily test.
2
The (t) curves should look di erent for large initial d =dt than for small initial d =dt. Why is this? 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 so these 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 nd that some of the trajectories form closed curves in phase space, while other trajectories do not. How does this relate to your results from part 2b?
Optional (but not di cult) 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.adding a 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 = np.linspace(0,-5,200) and pass this to odeint.
3.Newman Ex. 8.7: Cannonball with air resistance. The code you need is almost the same as for the Orbit of a planet” in-
