Computational Physics Lab 6

Initial Value Problems

This goal of this lab is to explore numerical modeling of the wave motion and to perform an exercise in stability analysis.

Part 1. Simple Wave Motion

Please download the file lab6_code.tar.gz and untar it:

tar xzvf lab6_code.tar.gz

This will create the directory "Lab6_code" and will place standard directory structure and example code files there. Compile the code for the initial value problem solving examples:

cd Lab6_code/src
make

Note that, if you are not running this on spudhammer, you might have to modify the file named "Makefile". You will need to find out where "numpy" is installed on your machine. On a Linux box this can be typically done as follows:

locate ndarrayobject.h | grep numpy | grep site-packages

Define the NUMPY_INC variable in the Makefile appropriately, using the directory name where you have found "ndarrayobject.h".

When you are done with the compilation, go to the ../tests directory and run the examples:

cd ../tests
setenv PYTHONPATH ../src
./LaxAdvectionSolver1d_example.py
./WaveSolver2d_example.py

The script "LaxAdvectionSolver1d_example.py" illustrates the motion of the advection wave using the Lax updating method discussed during the initial value problem lecture. If you find that the animation is too slow, try changing the matplotlib graphics backend to GTK:

./LaxAdvectionSolver1d_example.py -dGTK
./WaveSolver2d_example.py -dGTK

Change the time step of the advection equation solver: set the dt_factor variable on line 30 of file "LaxAdvectionSolver1d_example.py" to some other values and rerun the simulation. The value 1.0 in the example corresponds to the amplitude amplification factor of 1. Interesting values to use are 1.01 and 0.9. Can you explain what you see using the von Neumann stability analysis results for the Lax scheme (one of the lecture slides describes the analysis)?

Write you own C++ class for solving the wave equation in 1d (mechanical wave on a string). Look at the "LaxAdvectionSolver1d" example code and write something similar using the leapfrog updating scheme (formula 6.6 on page 158 of Giordano & Nakanishi). Please also look at the header file of the "CPInitialValueSolver" base class. Extensive comments in the header file explain the class usage. Note that you are switching from a 2 time steps scheme to a 3 time steps scheme, so the call to the base class constructor should be modified accordingly. Include support for three different types of boundary conditions: the ends are immobilized, the ends are free to move, and the ends are connected together (periodic boundary condition).

To build your code, modify the files "Makefile" and "CPInitialValue.i". The modifications are rather simple: if, let say, your class is named "WaveSolver1d" then insert the name of the object file "WaveSolver1d.o" into the Makefile together with "LaxAdvectionSolver1d.o", and use the name of the header file "WaveSolver1d.hh" in the SWIG wrap file "CPInitialValue.i" in the same manner as "LaxAdvectionSolver1d.hh" is already used.

Write a Python script which animates your wave equation solution (modify the "LaxAdvectionSolver1d_example.py"). Don't forget that you need to specify the initial conditions for two time steps (advection equation is first order in time, so it needs only one).

Explore the string behavior using different boundary conditions. You will see that the wave packets reflect from the boundaries, and they either change or do not change the phase depending on whether the ends are fixed or allowed to move. Can you explain why the advection waves do not reflect from the boundaries?

Explore different time steps. At what values of Δt the solution becomes unstable? Note how the shape of, let say, a triangular wave packet deteriorates with time in case Δt is less than |Δx/c| (for example, Δt = |Δ x/(2 c)|). Periodic boundary condition is the easiest to use for this study. Even though the von Neumann stability analysis shows that the amplitude amplification factor is 1, the system accumulates the phase error. Effectively, due to discretization errors wave components with different wavelength move with slightly different speeds. This phenomenon is called numerical dispersion -- you can find more details here or in the Numerical Recipes book.

Part 2. Stability Analysis

Apply the von Neumann stability analysis to the Lax-Wendroff differencing scheme for solving the advection equation:

This scheme is more precise than the Lax scheme: precision of each time step is O(Δt³), and it performs much better with small steps. Derive the stability criterion for this scheme.

Submit your lab report by email before 2 pm 03/04/2008. Provide the complete code necessary to solve the problems posed in Part 1 together with some plots of the results. Include a brief document (preferably in pdf format) with the derivation of the stability condition for Part 2.