780.20: 1094 Session 15

Handouts: "Three-Dimensional Plots with Gnuplot", excerpts from Landau/Paez PDE chapters, program printouts

Here we'll look at some basic examples of solving partial differential equations (PDE's) based on simple difference methods. Note: The PDE codes were designed to have very short printouts; there are many ways to improve the coding technique!

SET's: Course Evaluations

Thank you in advance for your comments, which are invaluable in planning future versions of this computational physics course. Some things you might comment on:

1. Please comment on the pace: too fast, too slow, not enough depth?
2. Do the background notes work? (Do you need them further in advance?)
3. What topics would you omit? Any ones you would spend more time on? Any topics not covered you think should be covered?
4. Do you think I should choose a text to go along with the course, even if we don't follow it closely?
5. Would it be better to spend more time on C++ programming (e.g., C++ classes)? If so, what would you cut to make time for this?
6. If you had a choice between computational physics courses taught entirely with Mathematica or MATLAB, entirely with C++, or entirely with Python, which would you choose? Or do you think a mix is more useful?
7. As always, suggestions for improvements to any aspect of the course would be appreciated.

3-D Plots with Gnuplot

We have frequently used Gnuplot for visualization, but we have only considered two-dimensional plots. Now we will want to make three-dimensional surface plots of functions and data (in particular, the output of today's programs!).

1. Follow through the handout on "Three-Dimensional Plots with Gnuplot". Make sure to look at the code for 3d_shape.cpp and predict what shape will be plotted.
   
   
2. Note that you can rotate a plot by grabbing it with the mouse.
3. Figure out how to make a parametric plot of a sphere using trigonometric functions. To plot a 2-d circle, you would type:

     gnuplot> set parametric
     gnuplot> plot [0:2*pi] sin(t),cos(t)

   (hint: think spherical coordinates). What command did you use for the sphere?
   
   
   

Laplace's Equation

Here we take a look at a solution to Laplace's equation by a simple relaxation method.

1. Take a look at laplace.cpp and see what problem it is solving and how it carries out the method described in section 25.4 of the handout. What are the basic approximations in implementing the PDE?
   
   
   
2. Compile and run the code and make a 3-d plot of the output using gnuplot. Include contours (try "set contour base" and "set contour surface"). To control the contour spacing, use "set cntrparam" to get contours at every 10 units (find out the details from "help set cntrparam"). What command did you use?
   
   
   
3. Given the equipotential lines, where is the electric field strongest?
   
   
   You might want to sketch in the electric field lines on your plot, starting from the 100V end.
4. How do you know if the relaxation method is converged? Pick a fixed point (i.e., fixed values of i and j) and see how much potl[i][j] changes from one iteration to the next (i.e., print out the value). Increase the number of iterations until you have five significant figure accuracy for your point. How can you build this into the code so that it stops iterating when this accuracy is reached? Describe a way to test the convergence more generally (you don't need to implement it).
   
   
   
5. Does it help to "damp" the iterations? Try adding damping (that is, the new value potl[i][j] at each (i,j) point is (1.-fraction) of the old value plus fraction times the formula for potl[i][j]). In fact, try 0 < fraction < 2 (which will reveal that "damping" is a poor choice of words!).
   
   
   
6. [Bonus] The numerical solution to the PDE can be used for any boundary conditions. Suppose you have designed a piece of equipment that is essentially a small metal box at 100 V within a larger, grounded one. You find that sparking occurs inside it, which indicates too large an electric field. You need to determine where the field is greatest so you can change the geometry and eliminate the sparking. Modify laplace.cpp (call it laplace_boxes.cpp) to satisfy the new boundary conditions and figure out where the electric field is most intense. How would you redesign the equipment to reduce the field?
   
   

Diffusion Equation

Here we look at a solution to the temperature diffusion equation that steps through in time using a (crude) finite difference method.

1. Take a look at eqheat.cpp and see how it carries out the method described in section 26.4 of the handout. [Note: We assume a cross section of 1 cm^2.] What are the basic approximations in implementing the PDE?
   
2. Compile and run the code to generate eqheat.dat. Look at eqheat.dat and then make a 3-d plot of it using gnuplot, using the comments in the code and the handout as guides. Include contours (which are isotherms here). Interpret the plot for your partner.
3. Stability test: Add a statement to the code to print out the constant so that you can compare it to 1/2 [note: the handout says 1/4 in (26.25) but this is a typo]. Then increase the conductivity to 0.4 and then 0.5. Look at the output file eqheat.dat before trying to plot (you should always do this!). Can you verify the stability criterion?
   
   
   
4. Now put the conductivity back to 0.12. Instead, let's change delta-x instead. How do we do that? What are delta-x and delta-t in the code? For what delta-x does the calculation become unstable?
   
   
   
   
5. [Bonus] Suppose instead of one 100 cm bar we have two identical bars, each 50 cm long, placed end to end. One bar starts at 100 degrees C while the other is at 50 degrees C. The outside ends are again put in an ice bath at 0 degrees C. Modify the code for this case and determine how the temperature varies with time and location with an appropriate plot.
   

Wave Equation

Here we look at the solution to the wave equation for a string with fixed endpoints.

1. Take a look at eqstring.cpp and see how it carries out the time-stepping method described in section 27.3 of the handout. Note how the wave speed and delta-x and delta-t are incorporated somewhat indirectly. What are the basic approximations in implementing the PDE?
   
   
   
2. Compile and run the code and make a 3-d plot of the output using gnuplot (some hints for gnuplot settings are in the comments). Estimate the wave velocity c from the graph (be careful of the scale of the time axis) and compare to the theoretically expected value.
   
   
   
3. The stability condition is said to be that the wave speed c is less than or equal to c' = Delta-x/Delta-t. Verify that the result is unstable if this is violated. What did you do and what happened?
   
   
   
4. Modify the code so that the initial conditions are that the string is plucked in the middle (rather than near one end). What did you have to change?
   
   
   How would you modify the code so that it is easier to implement more general initial conditions?
   
   
   

Bound States in Momentum Space: Delta-Shell Potential

As discussed in the background notes, the bound-state Schrodinger equation in momentum space can be solved as a matrix eigenvalue problem. We set it up so that the matrix is symmetric so we can either use the GSL routines we've used before. We'll also dynamically allocate some arrays using "new" and "delete".

The delta-shell potential is a convenient choice for a test potential because we can find the eigenvalues analytically and it takes a simple form in momentum space. It has only one bound state; if we only print out the negative-energy (E<0) eigenvalues then we should find just one. See the background notes for additional details. [Warning: It is not a well-behaved potential in momentum space, so you will be challenged to find good solutions.]

1. What are the units of lambda in the delta-shell potential? [Note: hbar=1 here.]
   
   
2. For mu=0.5, b=2.0, and lambda=-1.0, use Mathematica to find the correct eigenvalue for the (single) bound state in the delta-shell potential by using the appropriate equation in the background notes. [Hint: FindRoot] What is the answer?
   
   
3. Look at deltashell_boundstates.cpp. Note how the non-gsl arrays are dynamically allocated and removed using "new" and "delete". If we had an int called size and wanted to create an array called my_array with dimension size, what command would we give?
   
   
4. Compile and link deltashell_boundstates.cpp using the makefile. Run with lambda=-1.0 but specify on your own the number of gauss points to use and the upper bound for the integrations of momentum. The latter may be much larger than you would guess because the delta-shell potential is actually rather poorly behaved in momentum space! Based on how the potential falls off with k, what do you think is a reasonable momentum cutoff? (What should you do to discover how the potential falls off?)
   
   
5. Run the code with different numbers of gauss points and different momentum cutoffs to try to improve the predicted eigenvalue. Estimate the precision of your best value for the bound-state energy. How could you improve the numerical accuracy of the calculation? (For example, would it help to change "job" when calling "gauss"? [Hint, hint.])
   
   
   
6. What does the momentum-space wave-function look like (it is printed out to a file)?
   • Does it fall off at large k? Does it oscillate? Is it well behaved at the origin? Can you account for its characteristics by some simple arguments?
     
     
     
   • The analytic wave function in coordinate space is given in the background notes. Try to calculate the momentum space wave function from this (as discussed in the notes) using Mathematica.
   • Is the momentum wave function normalized? What is the normalization condition?
     
     
7. Bonus: Verify that, regardless of the potential's strength, there is only a single bound state.
8. Bonus: Modify the code to calculate the l=1 bound state.