*Handouts: "Three-Dimensional Plots with Gnuplot",
excerpts from Landau/Paez PDE chapters, program printouts,
"Using the GDB Debugger"
*

In this session we'll look at a some basic examples of solving partial
differential equations (PDE's) based on
simple difference methods and take a quick look at one
aspect of optimization.
**Note: The PDE codes were designed to have very short
printouts; there are many ways to improve the coding
technique!**

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!).

- 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.*

- Note that you can rotate a plot by grabbing it with the mouse.
- 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?*

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

- 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?*

- 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?*

- 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. - 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).

- 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!).

- [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?

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

- 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?*

- 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.
- 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?

- 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?

- [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.

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

- 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?*

- 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.*

- 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?*

- 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?

One of the most common floating-point operations is to square a number. Two ways to square x are: pow(x,2) and x*x. Let's test how efficient they are.

- Look at the printout for the square_test.cpp code. It implements these two ways of squaring a number. The "clock" function from time.h is used to find the elapsed time. Each operation is executed a large number of times (determined by "repeat") so that we get a reasonably accurate timing.
- We're going to set the optimization to its lowest value to
start. In Dev-C++ under Tools->Compiler Options->Compiler, add
`-O0`("minus oh zero") to the top box for the Compiler Command Line (if not already there). Check that Tools->Compiler Options->Settings->Optimization all say "no". - Compile square_test.cpp and run it.
Adjust "repeat" until the
minimum time for each is at least 0.1 seconds.
*Which way to square x is more efficient?*

- If you have an expression (rather than just x) to square,
coding (expression)*(expression) is awkward and hard to read.
Wouldn't it be better to call a function (e.g., squareit(expression)?
Add to square_test.cpp a function:

`double squareit (double x)`

that returns x*x. Add a section to the code that times how long this takes (just copy one of the other timing sections and edit it appropriately).*How does it compare to the others? What is the "overhead" in calling a function (that is, how much extra time does it take)? When is the overhead worthwhile?*

- Another alternative, common from C programming:
use
`#define`to define a macro that squares a number. Add

`#define sqr(z) ((z)*(z))`

somewhere before the start of`main`. (The extra ()'s are safeguards against unexpected behavior;**always**include them!) Add a section to the code to time how long this macro takes; what do you find?

- One final alternative: add an "inline" function called square:

`inline double square (double x) { return (x*x); };`

that is a function prototype**and**the function itself. Put it up top with the squareit prototype. Add a section to the code to time how long this function takes.*What is your conclusion about which of these methods to use?*

- Finally, we'll try the simplest way to optimize a code: let the
compiler do it for you! Change the compile flag
`-O0`(no optimization) to`-O2`(that's the uppercase letter O, not a zero). Recompile and run the code.*How do the times for each operation compare to the times before you optimized?*(If necessary, increase`repeat`to make all times above 0.05 seconds again.)

- In your project programs, once they are debugged and running,
you'll want to use the
`-O2`(or maybe`-O3`) optimization flag.

furnstahl.1@osu.edu