Exercise 3: Implicit Equations

This is an exercise in using IDL for program development, plotting and debugging (finding errors in) a program. These particular errors are "planted", but are typical of the types of errors you will almost certainly make when you write your own programs.

You will first get acquainted with basic use of IDL, then work mainly with IDL's development and debugging environment IDLDE.

The questions today are organised slightly differently. There are a number of links to "hints" that you may use if you get stuck. Often the answers were discussed during the lecture, and/or can be found by looking in the PDF lecture notes. However, the hints are not very expensive, so don't worry too much about using them.

Pick up your bonus points here!   Credits: 6/6

To update the exercise files (which includes getting the files for this weeks exercise), do

    cd ~/ComputerPhysics
    cvs update -d

In case of problems, see the CVS update help page.

The purpose of this part is to introduce you to some basic IDL commands and to IDLDE, the IDL Development Environment. Start the IDLDE, either from the dropdown menus, the icon on the desktop or by typing "idlde" from a command prompt. This opens a window with several frames. The only one relevant right now is the one titled Command Line, located at the bottom of the window. Use the on-line IDL Basics (This is an old version, but the basics have not changed since 1997) and work through at least the first Section (PDF pages 9-12, document page nos 3-6); trying out the examples on the screen

• While working on the examples, be sure to get used to using the keyboard arrow keys to get back to and modify old input lines.

• You may also want to try out the (UNIX/Linux) keyboard short cuts CTRL-a, CTRL-e, CTRL-b, CTRL-f, CTRL-k, CTRL-u, CTRL-z -- not to mention shift-CTRL-L that lists all available short cuts.

• Choose a string from a command that was executed some lines ago, type ^string<RETURN>, where ^ is the caret  ("hat") character. This takes you back to that particular string in the command history.

If you have time (within the 30 minutes), look at some of the other Sections in the manual as well. Some other time, browse the first few Chapters in Getting Started with IDL.

• "Browse" means "skim" / "look it through", reading only just enough to know what is there, so you know where to go back and read more carefully when the need arises. Browsing is a crucial technique in these Information Age days; there is 10^N times more information out there than you will ever have time to assimilate.

The roles of journal files (with names ending in .jou) and procedure files (.pro) are explained in the Lecture Notes, and in this hint:

Summary of IDL journal files and procedures   Credits: -1/-1

Here are some follow-up questions:

Let's assume 060915.jou is an IDL journal file, prog.pro is an IDL main program, doit.pro is an IDL procedure with two arguments, and fun.pro is an IDL function with three arguments. Indicate if the following are syntactically correct or incorrect, assuming the internal machinery knows how to handle it:      (scores +1 and -1)

IDL> @060915.jou correct incorrect Credits: 1/-1 OK

IDL> .run prog.pro correct incorrect Credits: 1/-1 OK

IDL> .r prog correct incorrect Credits: 1/-1 OK

IDL> @doit correct incorrect Credits: 1/-1 OK

IDL> doit,a,b correct incorrect Credits: 1/-1 OK

IDL> a=fun(a) correct incorrect Credits: 1/-1 OK

  From now on you will be working on material in the 3_Implicit directory, so you need to make sure that your IDLDE session is running in that directory. So,

• Be sure to set the option Window -> Preferences -> IDL -> editor -> Change directory on open (click the small triangle in front of "IDL" to see the "editor" options)!

• Open the file Implicit.pro in the ComputerPhysics/3_Implicit from the FILE menu or the folder icon in the toolbar and try to compile it. To compile, either
  • choose 'Compile Implicit.pro' from the 'Run' menu, or
  • press the 'compile' button in the toolbar -- the yellow gear wheels. There are tool tips (small pop-up texts) to help identify buttons

    

• You will discover that there are syntax errors (typical of what one often happens to introduce by mistake) in the file. IDLDE reports the line number of the syntax error in the Console window. Click on the link to get to the actual position in the editor window.
  
  • Fix the syntax errors, by editing the code in the IDLDE window panel.

• Please note: A "good" (but somewhat obsolete) version of the program is included with the Section 3 lecture notes, but it is rather pointless to immediately use it for answering the questions. It is not primarily the points of this Exercise that you need, but rather the training that comes with trying to find the errors. So, do yourself a favor, and use the program text in the lecture notes only as a last fall-back, if you get stuck for too long.

The Implicit.pro file contains one procedure, NewT, and two functions, Ne_over_N and residual:

NewT

computes factors that only depend on temperature

Ne_over_N

computes the relative number of electrons, given the electron density

residual

is the function that vanishes when a consistent solution has been found. The exact form of the function may be changed, as long as it returns zero for the same physical solution. By changing this function in various ways one can improve the speed (reduce the number of iterations needed).

We first test if the NewT procedure works as it should:

NewT problem 1

• Try calling the NewT procedure:

    IDL> newt, 3000.  

  The result will be a 'run time error' (unless the code has been corrected). The calculation of the variable c (see Eq. (9) in Section 3 of the Lecture Notes) has a problem that should be obvious from the run-time diagnostics.

  If that is not enough, here is a hint:

  Hint about undefined variable   Credits: -1/-1

NewT problem 2

• After fixing the problem with the undefined variable, try again

    IDL> newt, 3000.  
  

  the result is IDLDE diagnostics about 'Floating illegal operand'. To figure out where these occur, use the menu entry Window -> Preferences -> IDL -> Interpreter -> Report floating point math exceptions, setting it to "Immediately (2)" (the same results can be obtained by setting the IDL system variable !except=2, but that can be difficult to remember some other time). Retrying the command

   IDL> newt, 3000
  

  now should tell you the line number where the problem occurs -- but actually in the current version of IDLDE it doesn't.

  To take care of this for now, two lines containing check_math force the code to stop and enter debug mode on math errors.

  To also make the line numners visible in the editor, right-click in the left hand side margin of the editor window and choose "Show line numbers".

• The cause of this problem is discussed in the Section 3 lecture notes, so with the additional information you get about the line number you should be able to realize what the problem is.

  If knowing where the problem occurs is not enough, you may get more help by using a 'break point':

• Click in the gray margin to the left of the line numbers in the editor window at the line with c = ... in the NewT procedure, and then click on Toggle break point (also available in the "Run" menu and with key-board short-cut Shift-Ctrl-B). A small blue dot will appear

  To set or remove a break point you can also just double-click in the grey margin.

• Now try running again:

      IDL> newt, 3000
  

  A new window opens, asking if you want to change the layout to debugging mode. Answer Yes and the windows inside IDLDE will be rearranged. IDLDE stops at the break point in the NewT prodedure, and marks the place with a small blue arrow.

  At this point you have access to all the local variables at the command line, so you can try the calculation of c on the command line. Use 'copy and paste' from the code to the command line, looking at the expression piece-by-pice.

  The values of variables are also available by just "hovering over" (pointing at) variable names in the editor window. A small window should then appear, showing the value of the variable.

  You can also try executing the line by clicking the "Step over" button (yellow arrow) in the Debug window pane.

• If none of the above helps, as a last way out, take this hint (but you might as well save the credit and read the lecture notes):

  Hint about very small numbers   Credits: -1/-1

• Recompile the code after fixing it, and try again until there is no problem with that line!

NewT problem 3

• This should now work OK (disable the break point and try):

      IDL> newt, 3000
  

• Now try a lower temperature:

      IDL> newt, 2000
  

  IDLDE now probably reports (in the log panel)

  % Program caused arithmetic error: Floating overflow
  

  To find the cause for this we again use the break point: Enable the same break-point as before, and after the code stops there click the "Step over" button repeatedly. You will notice an overflow occuring at some point. One of the K_ variable is 'Inf' (infinity), which is the result returned when an arithmetic experession 'overflows' (becomes larger than about 3 10³⁸).

  Which of the K_ variables is it? K_H K_He K_Fe Credits: 1/-1 OK

  The cause and cure of this problem is also discussed in the lecture notes, so read about it, and then fix the code!

  In brief, what you need to do is this:

• Change the definitions of the all three K_Xx variables, so they become the inverse of what they were before
  • note that you need to rewrite the expressions in terms of negative exponentials, and not just invert the expressions after they have been computed. Remember to invert the prefactor as well!

• Having done that, you need to also change the expressions where the K_Xx variables are used, in the Ne_over_N function. There is a hint about how to do this, but this was also discussed during the lecture and is in the lecture notes so try without the hint.

  Hint about how to change the expression   Credits: -1/-1

    cd ~/ComputerPhysics/3_Implicit
    cvs diff

Ne_over_N tests

  Ne_over_N(N_e) = N_e/N

• Plot the function over the entire range, from N_e=1 to N_e=1e30, by entering

      IDL> newt, 3000.
      IDL> n_e = 10.^findgen(30)
      IDL> print, n_e
      IDL> plot, n_e, Ne_over_N(N_e), /xlog, /ylog
  

  A problem that should be ignored: This command causes an overflow in the IDL plot function (probably due to some minor stupidity in the calculation of tickmark values), which you should just ignore.

  NOTES:

  • IDL is case insensitive, so N_e and n_e are the same!
  • plot, .., /xlog, /ylog produces a log-log plot

  Over what range of function values does Ne_over_N() vary at this temperature (3000 K) with N_e defined as above? 10 0- 10-15 10 0- 10-30 10 0- 10-20 Credits: 3/-3 OK

• Add curves for other temperatures, by typing

      IDL> for T=4000,15000,1000 do begin $
      IDL>   newt,T & oplot,N_e,ne_over_n(N_e) & end
  

  NOTES:

  • the $ at the end of a line is a continuation character; you may type the whole sequence of commands on one line if the window is wide enough
  • the & (ampersand) allows more than one command to be typed on one line
  • the final end (matching the begin) is redundant in command mode

• Now add a set of lines corresponding to N_e/N, for various N from 10⁵ to 10²⁵:

      IDL> for lgn=5,25 do oplot, N_e,N_e/10.^lgn
  

  NOTE: no "begin ... end" is needed, since there is only one command after "do"

  So, on the screen you now have two families of curves, corresponding to various values for T and N. The solution for any given pair of values is given by the intersection of the respective curves.

• But the plot now looks pretty cluttered. Redo the last few commands (from plot,...), using the arrow keys, but load the rainbow color table, and colorize the curves, with

      IDL> device, decompose=0		; for 24-bit, direct color
      IDL> loadct, 39            ; rainbow colors
      IDL> for T=3000,15000,1000 do begin $
      IDL>   newt,T & oplot,N_e,ne_over_n(N_e),  color=!d.table_size*T/16000.
      IDL> for lgn=5,25 do oplot, N_e,N_e/10.^lgn, color=!d.table_size*lgn/26.
  

  This should give you colored families of curves, where you can more easily see which values are low (blue) and high (red).

  • Explanation: Depending on circumstances (e.g. how many colors the display station has, and how many that are already taken by other applications), the number of colors available to IDL may vary. The variable !d.table_size contains the number of colors in the table; here it is used to make color choices independent of the table size.

• From these curves you should already be able to read off solutions.

  What, for example is the closest N_e/N for T=3000 and N=1e22? 1e-1 1e-4 1e-6 Credits: 1/-1 OK

  To make it easier to identify the relevant curves, you may want to overplot the curves for T=3000 and N=1e22 in thick white, by doing

      IDL> newt, 3000
      IDL> oplot, N_e,ne_over_n(N_e), thick=2
      IDL> oplot, N_e,N_e/1e22, thick=2
  

• Before continuing, set the !except variable back to the default value (1), to avoid a lot of output about (non-critical) arithmetic underflows:

      IDL> !except=1
  

• You should now be able to run the main program, by clicking first the 'Compile' and then the 'Start execution' buttons.

  

  When the program runs, it prints the accumulated number of iteration in zbrent and makes a surface plot of the result.

  How many iterations per value does it use on the average (choose the closest value)? 5 10 15 20 Credits: 1/-1 OK

• As is, the program requires an unnecessarily large number of iterations. Change the code to work in terms of the logarithms of N and N_e.
  • But before you do, add three lines with plotting commands to the code. One,

        plot, [1,1e30],[1e-10,1], /xlog, /ylog, /nodata
    

    in the main program, just after the call to NewT, and two,

        plots, N_e, N_e/N, psym=1, symsize=0.5, col=0.6*!d.table_size
        plots, N_e, Ne_over_N(N_e), psym=1, symsize=0.5, col=0.75*!d.table_size
    

    inside the residue() function. "plots" plots a symbol (1 is "+"), every time the residue function is called while the zbrent procedure is looking for a solution to N_e/N = Ne_over_N(N_e), and through the color coding you can distinguish which family of curves the points belong to (cf. the earlier plots).

  • A pattern of convergence flashes by on the screen as you run the code, and through it, you should be able to appreciate the effect of the changes introduced by going over to logarithmic variables. However, the plots are probably flashing by too fast, so you should add

        wait,0.02		; wait 0.02 seconds
    

    or something similar, just after the "plots" calls; a 0.02 second wait should actually be enough to see points being plotted sequentially

• Do the change to logarithmic variables in two steps; first changing the return value in the "residue" function from 1-f to alog(f). This corresponds to changing to a logarithmic y-axis in the f(x)=y problem.

  How many iterations does it take on the average now? 6 9 14 Credits: 3/-3 OK

• Then change the "x-axis" as well, by changing the name of the argument to residue from "N_e" to "lg_N_e"; i.e., you now assume that the argument is the logarithm of "N_e" (base-e or base-10 after your own taste). This is a bit more tricky than changing the y-axis, so consider the points below carefully:
  • "N_e" is still to be used in the expression, so you need to add a line that calculates it from "lg_N_e"
    • Note that, unfortunately "ne" is not a valid variable name.  NE is one of the reserved words used by IDL (means 'not equal' in a logical statement) and, since IDL is case insensitive, ne is reserved as well.

  • If "residue" is supposed to have a logarithmic argument, then the starting values given to "zbrent1" must be logarithms as well,

    OK? yes no Credits: 1/-1 OK

    Change this directly in the argument list, at the line where zbrent1 is called.If you don't understand how, take a look in the lecture notes, or use this hint:

    Hint about how to change the argument list   Credits: -1/-1

  • Also, the return value of zbrent1 is "that x-value for which residue(x) vanishes", hence it too is a logarithm. But the result that goes into the array f_e should still be the same as before.

    Make one more change of the same line (the zbrent1 call) that takes care of this (there is another hint here, but you may also consult the lecture notes).

    Hint about another change   Credits: -1/-1

• If you have succeeded in changing variables correctly, you should notice a significant decrease in the number of iterations.

  How many iterations does it take on the average? 5 8 11 Credits: 3/-3 OK

• If you have time, think about ways of improving the initial values.
  • How about using the result from the previous Section, where only hydrogen was ionizing?

    If you adopted that value of N_e/N (from the previous exercise) ...

    ... would it be smaller larger Credits: 1/-1 OK

    than the answer? If you have time, implement it as one of the initial values.

    How many iterations now? 3 4 5 Credits: 3/-3 OK

  • Another, even simpler idea is to use the value of N_e obtained at one pressure as the initial guess for the next pressure. Since the pressures are increasing each time through the innermost loop,

    should you use the previous value for Ne_min Ne_max Credits: 3/-3 OK

  • Any idea for a better Ne_max value? This is a difficult question, and you may choose to ignore it, but if you want to give it a try, here is a free hint: Look at the slopes of the red and yellow families of curves (on the plot you made earlier -- just redo the previous commands if the plot is gone).

Spent a couple of hours of follow-up on the Lecture, studying the chapter in Numerical Recipes about root finding, and browsing the Fortran book.

After the exercise, it would be a good idea to continue reading some of the introductory sections of the Getting Started in IDL and IDL Basics manuals and the Building IDL Applications manual. Go through the material from the exercise; making sure you understand what is going on, by going back and forth between reading the manuals and the program text.

You may also want to try out some of the examples on screen.

If you're interested in history/stories; here is the story behind Murphy's law.