PHYCS 6730 Lab Exercise: Crank Nicholson Method

To integrate a parabolic partial differential equation, such as the Schroedinger equation or the diffusion equation the Crank Nicholson method provides stability and a result good to second order. For the diffusion equation 
(1)   dP(x,y,t)/dt = alpha2 (d2/dx2 + d2/dy2)P(x,y,t)
we can discretize in space and time so that the time derivative is written either as a forward difference (explicit) 
(2)   P(x,y,t+dt) = P(x,y,t) + alpha2 dt D P(x,y,t),
where D is the discrete Laplacian, or backward difference (implicit) 
(3)   P(x,y,t+dt) = P(x,y,t + alpha2 dt D P(x,y,t+dt)
Adding (2) and (3) and solving for P(x,y,t+dt) gives 
(4)   (1 - alpha2 dt D/2)P(x,y,t+dt) = (1 + alpha2 dt D/2)P(x,y,t)
This is the Crank-Nicholson form. Stepping forward in time involves solving a sparse linear system.

The code ~p6730/exercises/diffusion/diffusion.cc solves the diffusion equation in one dimension with boundaries at infinity. This problem is quite easily solved with a Fourier transform, but in coordinate space the matrix D for the second derivative generates a tridiagonal matrix, which is also very easy to solve using the standard Crout reduction.

Exercise 1. Understanding the code

The code that solves the problem is in diffusion.cc. Read through the code and answer these questions:
  1. What is the value of the function P(x,0) for all x (at time t = 0)?
  2. What is the purpose of introducing a new class PVec?
  3. From what class does the PVec class inherit its properties?
  4. What properties of the NRVec class are changed in PVec?
  5. Why is the LU decomposition done outside the loop over t?
  6. Write the definition for a class B1Vec derived from the base class NRVec that uses 1 as the base for subscripts, rather than the conventional 0.
Notice that we used the public functions in the NRVec class to access the vector component and to determine the length of the vector. A simpler approach would have been to say return v[x+nn/2], but the members v and nn are declared private in the NRVec class. So our derived class can't use them. If they had been declared protected instead of private, we could have. But if we then later changed the names of those members in the base class, our derived class would break.

Exercise 2. Running the code

Use the code to solve the same diffusion equation problem as in the homework assignment this week. That is, solve for alpha2 dt = 0.1 with 25 time steps. Check that the total probability is conserved. Plot and compare with a Gaussian. Estimate the standard deviation from the 2/3 rule. Is it consistent with what you expect from the analytical solution in the notes?