Physics and steepest descents: sch.C

 

Write a program called sch.C based upon the code fragment below (and your physics.C program from HW#9) to solve Schrödingers's equation for the three-dimension harmonic oscillator. Use the following parameters:

• n=16 (size of FFT grid)
• d=0.75 bohr (grid spacing)
• k=0.25 hartree/bohr (spring constant)

Hint: In quantum mechanics, the analytic solution for the energy of the three-dimensional harmonic oscillator is where . In our units , so the expected result is 3/4 hartree. With the above resolution (d) and box size (n), the numerical result which you should find in the end is . (The fact that d and E are equal is just a coincidence in this problem. You can try running with different d's and see how the results for E change.) After about 400 iterations of the steepest descent algorithm (with alpha=0.0004 as defined below), you should reach an energy within about 1 micro-hartree of this final result.

Note: The algebraic results derived/to be derived in lecture for the energy and gradient are

where

and are the sample values of the potential on the real-space grid points, .

Code fragment:

                        .
                        .
                        .
const int Itmx=400;
main() {
                        .
                        .
                        .
  // Open file to store convergence results
  FILE *fp;
  if ( (fp=fopen("conv.dat","w"))==NULL ) {
    printf("Unable to open file!\n");
    exit(1);
  }
                        .
                        .
                        .
  n=16; d=0.75;  // Parameters for the calculation
                        .
                        .
                        .
  V.evalr(sho,0.25); // SHO potential stored in pvector V

  Vtilde= ... // See formula above

  // Initialize with normalized, but otherwise random, wave function
  u.eval(rnd); // Random elements
  u=u*(1./sqrt(u^O(u))); // Normalize u for starting
                        .
                        .
                        .
  alpha=4e-4; // Step size for going down the gradient
  for (int i=0; i<Itmx; i++) {
    // Compute energy and gradient for wave function u
    Denom= ... // See formula for ``D'' above.  Hint: don't forget to use real()
    Num= ... // See formula for ``N'' above.  Hint: use definition of Vtilde
    E=Num/Denom;
    grad= ...

    // Output: 
    //   iteration number, energy, length of u, cosine of angle between grad and last dir
    printf("%4d %20.15f %e %e\n",
	   i,E,sqrt(Denom),real(grad^dir)/sqrt(real(grad^grad)*real(dir^dir))
	   );
    fprintf(fp,"%4d %20.15f %e %e\n",
	   i,E,sqrt(Denom),real(grad^dir)/sqrt(real(grad^grad)*real(dir^dir))
	   );
    fflush(fp);

    // Choice of search direction (opposite to the gradient, for now...)
    dir=(-1.)*grad;

    // Put code for computing best alpha between here...

    // ... and here when you do Part 3.2.


    u=u+alpha*dir; // Take a step of size alpha down along the search direction
    printf("         ===> alpha=%e\n",alpha); // Remind us of step size 

  }

}

• Gathering data