# 780.20: 1094 Session 7

Handouts: eigen_basis_class.cpp and diffeq_oscillations.cpp printouts.

Now that we've got routines to solve differential equations, we're going to explore some interesting ones: nonlinear oscillators. Today we'll play with a program that solves for the time dependence of such an oscillator.

Your goals for today (and ...):

• If you didn't complete it, do the plot from Session 6 of relative error at t=1 vs. mesh size h.
• Think about how to enhance the eigen_basis code with more C++ classes.
• Run a code that solves the differential equation for a (driven) nonlinear oscillator and explore how the time dependence changes as various input parameters change.
• Add friction (damping) to the code.

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

## (Possibly) Leftover Task from Session 6

Spend about 45 minutes (or less) on this.

1. Integrating a First-Order Differential Equation. Try to finish through part 7. If you find that the errors for Euler and Runge-Kutta lie on top of each other, most likely you have not evaluated the exact answer at precisely the same time as the last points. If you get stuck, ask an instructor.

## More on C++ Classes: eigen_basis_class

The code eigen_basis_class.cpp is a simple modification of eigen_basis.cpp to use the Hamiltonian class we introduced for eigen_tridiagonal_class.cpp. Here we'll take a few minutes to think about how to introduce additional classes.

1. Take a look at the eigen_basis_class.cpp printout and note how the Hamiltonian class is re-used without modification. (If you haven't done so yet, read the discussion of this class in the Session 7 notes.) The only tricky change is that matrix indices go from 1 to dimension rather than from 0 to dimension-1. What parts of the Hamiltonian class implementation do you not yet understand?

2. The potential is another good candidate for a class. We'd like to just evaluate the potential at r without having to use constructions like the switch statement in the Hij_integrand function with all the messy void parameters. (Think about how awkward and prone to error it is to add another potential.) What would you like the declaration statement for the Potential class to look like? What method(s) would you like the class to have?

3. Give at least one example of an additional class that would be useful to define.

## Driven Nonlinear Oscillations

The Session 7 notes describe the driven nonlinear oscillator that is coded in diffeq_oscillations.cpp. Note that the force depends on k and an exponent p, the external force has a magnitude f_ext, a frequency w_ext, and a phase phi_ext. The initial conditions in position and velocity are designated x0 and v0. You also have control over the time interval (increase t_end to see longer times), the step size h, and how often points are printed to the file (plot_skip).

1. Use make_diffeq_oscillations to create diffeq_oscillations. This code outputs to the file diffeq_oscillations.dat five columns of data: t, x(t), v(t), kinetic energy, and potential energy. There are four gnuplot plot files provided (diffeq_oscillations1.plt, etc.), each of which generates a different type of plot. Run diffeq_oscillations with the default values (enter "0" when it says "What do you want to change?") to calculate a data set. Start gnuplot and "load diffeq_oscillations1.plt" and then "load diffeq_oscillations2.plt". (Once you've given these commands once, you can use just use the arrows to go back and forth.) Briefly, what do each of these plots show?

2. Wouldn't it be convenient to generate all four plots at once in separate files? Load "diffeq_oscillations_all.plt"!
3. It's always a question whether or not you have coded a problem correctly, so you should always seek ways to check your results. One possibility is if we have a known solution. This works for p=2 (simple harmonic oscillator). What about other p? Another check is to identify a quantity that shouldn't change with time. Create a plot of such a quantity (you'll want to increase t_end) and observe the effect of changing the step size h to a larger value [e.g., try 10 and 100 times larger]. How do you decide on a reasonable h to use?

(The "plot_skip" parameter indicates how often a point is written to the output file. So plot_skip=10 means that every 10 points is output.)
4. Verify that different amplitudes (e.g., different initial conditions determined by x0 and v0) lead to different periods for an anharmonic oscillator (p<2 or p>2). [Hint: You might find the "append" option useful.] Can you identify a qualitative rule? E.g., does larger amplitude mean shorter or longer period always? Can you explain the rule?

5. Go back to the original parameters (quit the program and start it again), which has p=2. Now add a driving force f_ext=10 with w_ext=1 and look at the time dependence and phase-space plots. Then increase w_ext to 3.14 and then to w_ext=6.28. What are you observing? Now repeat with p=3 (starting with f=0). Can you find resonant behavior?

Real-world systems have friction, which means the motion will be damped. The Session 7 notes have a list of three simple models for friction. We'll implement viscous damping: Ff = -b*v, where v(t) is the velocity.

1. Introduce the damping parameter "b" into the code:
1. add it to the force_parameters structure (with a comment!);
2. add it to the list of local force parameters in the main program;
3. give it an initial value;
4. add a menu item (e.g., [13]) and a case statement to get a new value.
Try this part out before proceeding.
2. Modify the "rhs" routine to include damping (you're on your own here!). What did you add?

3. Test your routine starting with p=2 and a small damping and look at both the time dependence and the phase-space plots. Then try some other p values.
4. Identify the three regimes described in the Session 7 notes: underdamped, critically damped, and overdamped.

## EXTRA: Using TWiki

By now, most everyone knows about wiki's, wikipedia, and the usefulness of having a web page that multiple users can interactively edit and post content to. Whether your research is theoretical or experimental, this can be an extremely useful tool for organizing a research project where you and a few other people need to look at plots and comment on them, or perhaps there is some experimental procedure that can be documented with a wiki and edited on the fly by multiple folks. To encourage this kind of online note-keeping, the folks in the Physics Computing Facility have set up a wiki capability. Here we'll try it out.

1. Go to http://www.physics.ohio-state.edu/TWiki/bin/view/Physics780/WebHome and register using the link in the upper left. Did you have any problems?

2. When you've successfully registered, return to the Physics 780 TWiki page and look at the ChrisProject example project. Click the edit link (on the left of the bottom menu) to see the format.
3. Now return to the Physics 780 TWiki page and try to create your own page (which you can use for your project!). Do this by clicking "Edit" on the main Physics 780 Twiki page. You'll see a section like this:

---++ %MAKETEXT{"Project Pages"}%
* ChrisProject Example Project page
* ...
* ...

By one of the asterisks put in the name of your page, e.g., LastnameProject or FirstnameProject. TWiki is set up to recognize text with this format (a single phrase with two caps and no spaces) as a link to a new page with that title. Save and see if you see something like this on the main site:

LastnameProject?

Click the question mark (in TWiki of course!) and a new page should start with that title. Notice that you can attach files, by clicking "Attach" (the file attach.txt is there as an example). These can be gnuplot plot files or figures or whatever. Note that there is a way to attach plots (e.g. png or jpeg format) and have them appear in the page. Unfortunately, it's a little more complicated than it needs to be. Ask Chris about how to do this. Finally, notice that the page has a revision history. You can see earlier versions of the example project page by clicking on the History links in the bottom menu.

Feel free to spin off multiple pages for your project if it's useful, and if you think you could make use of something like this for your (summer) research the folks in the physics computing facility would be glad to create an entirely new page. Contact Bryan Dunlap in the physics computing facility (ground floor PRB) if you think this would be something you or your research group would find helpful.

## EXTRA: Looking for Chaos (Part I)

Now we want to put it all together: a damped, driven, nonlinear oscillator. A different system with the same basic features is the realistic pendulum, which is described in the Session 7 notes.

1. In the notes there is a list of characteristic structures that can be found in phase space, with sample pictures. Can you find combinations of parameters that produce pictures like these? (Try to imitate the x(t) vs. t pictures first.)