Theory:

1. With gnuplot/pyxplot how do you plot different sets of data (with different styles and titles) from the same data-file?
   Hints:
   1. manual → commands → plot → data → index;
   2. in gnuplot interactive mode: help index;
   3. in pyxplot interactive mode: help data index;
2. Suppose an array is declared as
   double a[] = {1,2,3};
   • What is the type and the value of a ?
   • What is the type and the value of *a ?
   • What is the type and the value of a+1 ?
   • What is the type and the value of *(a+1) ?
   • What is the type and the value of *(1+a) ?
   • What is the type and the value of 1[a] ?
   • What is the type and the value of a[1] ?
3. Given double x=5; what is the type and the value of the expression *(&x) ?

Practice: Quantum particle in a box or Standing wave in a string.

• Introduction:

  The wavefunction ψ(ξ) of a particle in a box satisfies the Schrödinger equation

  -ψ''=εψ ,

  where the prime denotes the derivative over the dimensionless variable ξ=x/L, x is the coordinate of the particle and L is the size of the box, ε is the dimensionless energy. The actual energy is E=(ℏ²ε)/(2mL²).

  The boundary condition is ψ(0)=ψ(1)=0.

  Numerically one can represent the wave-function on a grid

  ξ_i = ih, h=1/N, i=0...N,

  where the wave-function is represented by a vector

  ψ_i = ψ(x_i), i=0...N.

  The second derivative operation can then be approximated using the three point finite difference formula,

  ψ''_i = (ψ_i-1-2ψ_i+ψ_i+1)/h² .

  Now the Schrödinger equation becomes a matrix eigenvalue equation

  ∑_j=1,N-1 H_ij ψ_j = ε ψ_i, i=1...N-1,

  where—since the wave-function at the end-points is zero—the summation goes over 1...N-1 and vanishing components ψ₀ and ψ_N have been omitted.

  The matrix H in the three-point approximation is given as (taking into account the zero boundary condition and omitting the first and the last component of the wave-function vector)

  	(	-2	1	0	0	0	)
  	(	1	-2	1	0	0	)
  	(	0	1	-2	1	0	)
  -h²H =	(	...	...	...	...	...	)
  	(	0	0	1	-2	1	)
  	(	0	0	0	1	-2	)

  Solving the matrix eigenvalue equation amounts to matrix diagonalization.

• Exercise:
  1. Build the Hamiltonian matrix h²H as a gsl_matrix structure. Hint:

     int n=20;
     gsl_matrix *H = gsl_matrix_calloc(n,n);
     for(int i=0;i<n-1;i++){
       gsl_matrix_set(H,i,i,-2);
       gsl_matrix_set(H,i,i+1,1);
       gsl_matrix_set(H,i+1,i,1);
       }
     gsl_matrix_set(H,n-1,n-1,-2);
     gsl_matrix_scale(H,-1);
     
     

  2. Diagonalize the matrix using GSL subroutines for diagonalization of real symmetric matrices. Hint:

     gsl_eigen_symmv_workspace* w =  gsl_eigen_symmv_alloc (n);
     gsl_vector* eval = gsl_vector_alloc(n);
     gsl_matrix* evec = gsl_matrix_calloc(n,n);
     gsl_eigen_symmv(H,eval,evec,w);
     

  3. Sort the eigenvalues and eigenvectors.

     gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_VAL_ASC);
     

  4. Check that your energies are correct. Hint: the ratio should be close to unity (why?):

       for (int k=0; k < n/3; k++){ /* print first n/3 eigenvalues (rescaled) */
         double exact = pi*pi*(k+1)*(k+1)/(n+1)/(n+1);
         double ratio = gsl_vector_get(eval,k)/exact;
         fprintf (stderr, "n = %i E_n/E_exact = %g\n", k, ratio);
       }
     

  5. Plot the several lowest eigenfunction (and compare with the analytical result).

     for(int k=0;k<3;k++){
       printf("%i %f\n",0,0.0);
       for(int i=0;i<n;i++) printf("%i %f\n",i+1,gsl_matrix_get(evec,i,k));
       printf("%i %f\n",n+1,0.0);
       printf("\n\n");
       }
     

  6. You may try to investigate how the accuracy changes with the number of grid points.
  7. You may try impove the accuracy of the method by using the five-point formula for the finite difference approximation to second derivative from here.