Applied Physics 461b/861b

COMPUTATIONAL PHYSICS

Assignment 2

due Wednesday, January 29, 1997

Practical information

  1. Use of functions and subroutines: read p.17-18 and p. 679-680 in GT.
  2. Compiling and linking subroutines in files other than the main program: You can compile an individual subroutine using the command f77 -c subroutine.f. This creates an object file called subroutine.o. When you compile the main program and create the executable, you don't have to recompile each subroutine. Instead, use the command f77 -o executablename main.f subroutine1.o subroutine2.o. If appropriate, you can also compile everything at once using the command f77 -o executablename main.f subroutine1.f subroutine2.f.
  3. How to include a subroutine from Numerical Recipes: The source codes for the various subroutines are available on mercury in the directory /secf/src/recipes_f/recipes. When compiling, refer to this file by its full name, i.e. f77 -c /secf/src/recipes_f/recipes/subroutine.f. This will create an object file in the current directory, which you can refer to in linking simply by subroutine.o.
  4. Plotting subroutines for fortran programs: There is a package called PGPLOT available on mercury which allows you to create plots as part of your fortran program, thus eliminating the need to write out the plotting data and read it into some other program like Mathematica or CricketGraph. You should be able to preview your plots by using an xwindows emulator (Exodus on the Macintoshes), and can print them out by writing them to a postscript file and printing to a networked printer. It is easiest to learn by being giving a simple program such as simple.f that you can compile, run, and eventually modify for your particular plot. Compiling and running requires some tedious setting of environment variables and library calls; all the necessary commands are included in the file pgplot.sh which you can run by typing csh pgplot.sh filename where filename is the name of your fortran program file (e.g. simple). A handout describing the commands in this program is also available.

1. Order and accuracy of the Euler algorithm

(a) Use the program euler.f to solve the differential equation dy/dx=-xy with boundary condition y(0)=1 (this is another equation that can be solved analytically). Make a table of y(1) and y(3) as a function of h. Explain how the results confirm that the Euler method has a error of O(h).

(b) ``A simple and often stringent test of an accurate numerical integration is to use the final value of y obtained as the initial condition to integrate backward from the final value of x to the starting point. The extent to which the resulting value of y differs from the original initial condition is then a measure of the inaccuracy'' (Koonin, p.28). Use Mathematica to plot the error in y(0) on forward and backward integration to x=1 (i) as a function of h; and (ii) as a function of the actual error in y(1). Use these plots to comment on how you would use forward and backward integration to test the accuracy in a case where you do not know the exact answer.

(c) Modify the program to replace the Euler update with a call to the Runge-Kutta subroutine rk4 from Numerical Recipes. Repeat parts (a) and (b) for the new algorithm.

2. Coffee cooling

Translate the program cool on p.22 into fortran (for your convenience, the True Basic file is available here). Analyze the coffee cooling problem using problems 2.1 and 2.2 as a guide. You need not follow every instruction literally. For example, in 2.1.a, you may want to check your program by comparing to the analytical answer (being aware that the agreement won't be absolutely exact) and only resort to the ``hand computer'' approach if you have a bug you are trying to track down or are very unsure of your coding. As another example, in 2.1.c, you might want to try a computer plot (pgplot, mathematica, or whatever seems most efficient) instead of a hand plot. In writing up your work, give some thought to the clearest and most useful format (see Appendix 1A at the end of Chapter 1 for ideas).