780.20: 1094 Session 9

Handouts: Printout of nonlinear.nb notebook, excerpt from Landau/Paez Chapter 5, printout of spline_function2 (.cpp, .h, _test.cpp), "Using GSL Interpolation Functions", ode_test.cpp printout

In this session, we'll look at using a GSL adaptive differential equation solver and then look at interpolation.

Your goals for today and ...

• Finish up left-over Session 8 tasks, plus look at some additions.
• Try out a GSL differential equation solver.
• Try out various GSL interpolation functions on a simple example.
• Represent a sample wave function using cubic splines and test how many points are needed.

Please work in pairs (more or less). The instructors will bounce around 1094 and answer questions.

Follow-up to Session 8

Work on these tasks for the first 30 minutes of the session only, then move on.

1. Complete the Session 8 sections on "Damped, Driven Pendulum" and "Looking for Chaos".
2. The Mathematica notebook nonlinear.nb looks at the same type of analysis as diffeq_pendulum.cpp only it uses the Duffing equation. Look through it and follow along with the printout. This exercise is mostly just to expose you to how to do lots of useful things with Mathematica, for future reference. Any questions about Poincare sections or the power spectrum?
   
   
   
3. If you're familiar with Mathematica and have time, you can convert nonlinear.nb to study the pendulum instead. For now, you can directly run the "answer" notebook called pendulum.nb.
   
   

GSL Differential Equation Solver

The program ode_test.cpp demonstrates the GSL adaptive differential equation solver by solving the Van der Pol oscillator, another nonlinear differential equation (see the Session 9 background notes for the equation).

1. Take a look at the code and figure out where the values of mu and the initial conditions are set. Change mu to 2 and the initial conditions to x0=1.0 and v0=0.0 (y[0] and y[1]). Note the different choices for "stepping algorithms", how the function is set up and that a Jacobian is defined, and how the equation is stepped along in time.
2. Create a project ode_test and add ode_test.cpp. Compile and run. Create three output files using the initial conditions [x0=1.0, v0=0.0], [x0=0.1, v0=0.0], and [x0=-1.5, v0=2.0] (just change values and recompile). Notice how we've used a stringstream to uniquely name each file.
3. Use gnuplot to make phase-space plots of all three cases on a single plot. Print it out and attach it. What do you observe? This is called an isolated attractor.
   
   
   
   
   

GSL Interpolation Routines

We'll use the Landau/Paez text's example of a theoretical scattering cross section as a function of energy to try out the GSL interpolation routines. The (x,y) data, with x-->E and y-->sigma_th, is given in Table 5.1 (see the Chapter 5 excerpt).

1. Create a project and load the gsl_cubic_spline_test.cpp code. Instead of the sample function, you will change the program to interpolate the Table 5.1 data from the handout. Set npts and the (x,y) arrays equal to the appropriate values when you declare them. Declare them on separate lines. An array x[4] can be initialized with the values 1., 2., 3., and 4. with the declaration:
   double x[4] = {1., 2., 3., 4. };
2. Modify the code to generate a cubic spline interpolation for the cross section from 0 to 200 MeV in steps of 5 MeV. Output this data to a file for plotting with gnuplot.
3. Now generate the corresponding linear and polynomial interpolations. You can generate these all at once (e.g., by making three splines in succession) or by running the code three times and changing the interpolation type between compilations.
4. Generate a graph with all three interpolations plotted. For each curve, how does the resonance energy E_r (the peak position) and gamma (full width at half maximum) as seen on the graph compare to the theoretical predictions of (78,55)?
   
   
   
   

Cubic Splining [finish for PS#4]

Here we'll look at how to use cubic splines to define a function (in the form that GSL routines use) from arrays of x and y values. A question that always arises is: How many points do we need? Or, what may be more relevant, how accurate will our function (or its derivatives) be for a given spacing of x points?

1. Take a look at the printouts for spline_function2. In spline_function2.cpp, the GSL calls have been packaged into functions. Notice that spline_funct has the same arguments (double x void *params_ptr) as a GSL function. Can you see how spline_funct works?
2. The test program splines two different sets of arrays (with two different spline structures). Note how params_ptr is set equal to the address of the relevant spline_struct (which is declared in spline_function2.h). Make a project and run the program (it's not very exciting!).
3. Modify the code so that it splines the ground-state hydrogen wave function: u(r) = 2*r*exp(-r) (You can delete the second spline.)
4. Your job is to determine how many (equally spaced) points to use to represent the wave function. Suppose you need the derivative of the wave function to be accurate to 0.1% for 0 < r < 4 (absolute, not relative error) Devise (and carry out!) a plan that will tell you the spacing and the number of points needed to reach this goals. What did you do?
   
   
   
   
5. Now suppose you need integrals over the wave function to be accurate to 0.01%. Devise (and carry out!) a plan that will tell you the spacing and the number of points needed to reach this goals. (To try out integrals, use one of the GSL integration routines on an integral involving u(r) that you know the answer to; there's an obvious one!)