| [ | 6.13 | -2.90 | 5.86 | ] | [x0] | [6.23] | |
| [ | 8.08 | -6.31 | -3.89 | ] | [x1] | = | [5.37] |
| [ | -4.36 | 1.00 | 0.19 | ] | [x2] | [2.29] |
Solve it using GSL.
Hint: you can solve linear systems using one of the many methods, for example,
gsl_linalg_HH_solve
gsl_linalg_LU_decomp gsl_linalg_LU_solve
gsl_linalg_QR_decomp gsl_linalg_QR_solveSee also:
Consider a quantum particle in a box (or a standing wave in a string).
Introduction (you might also want to read the articles "Quantum particle in a box" and/or "Standing wave in a string"):
The wave-function \(\psi(x)\) of a particle in a box of length \(L\) satisfies the Schrödinger equation \[ \hat H\psi = E\psi \;, \] where \(E\) is the energy of the particle and where the Hamiltonian \(\hat H\) is given as \[ \hat H=-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} \;, \] where \(m\) is the mass of the particle. The boundary condition for an infinitely deep box is that the wave-function vanishes at the borders, \[ \psi(0)=\psi(L)=0 \;. \]
Let us introduce dimensionless variables \(\xi\doteq \frac{x}{L}\), \(\epsilon \doteq \frac{2mL^2}{\hbar^2}E\), and the new function \(u(\xi)\doteq\psi(\xi L)\). The Schrödinger equation then reduces to \[ \hat H' u = \epsilon u \;,\] where \(\hat H' = -\frac{\partial^2}{\partial\xi^2}\), with the boundary condition \[ u(0)=u(1)=0 \;. \]
Numerically one can represent this problem on an equidistant \(N\)-point grid with step \(s\), \[ \xi_i = is \,,\; s=\frac{1}{N} \,,\; i=0\dots N \;,\] where the wave-function is represented by a vector \(u_i = u(x_i)\), and where the second derivative operation in the Hamiltonian can then be approximated using the three point finite difference formula, \[ \frac{\partial^2u_i}{\partial\xi^2} = \frac{u_{i-1}-2u_i+u_{i+1}}{s^2} \;. \] Now the Schrödinger equation becomes a matrix eigenvalue equation, \[ \sum_{j=1}^{N-1} H'_{ij} u_j = \epsilon u_i \,,\; i=1...N-1 \;, \] or, in matrix notation, \[ H'u=\epsilon u \;, \] where the matrix \(H'\) is given as \[ H' = -\frac{1}{s^2}\left( \begin{array}{cccccc} -2 & 1 & 0 & 0 &\dots & 0 \\ 1 & -2 & 1 & 0 &\dots & 0 \\ 0 & 1 & -2 & 1 &\dots & 0 \\ \vdots &\vdots &\vdots &\vdots &\vdots & \vdots \\ 0 &\dots & 0 & 1 & -2 & 1 \\ 0 &\dots & 0 & 0 & 1 & -2 \\ \end{array} \right) \;. \] Notice that—since the wave-function at the end-points is zero—the summation goes over \(1\dots(N-1)\) and the vanishing components \(u_0=0\) and \(u_N=0\) have been omitted and the boundary condition is thus built in the formulation of the matrix problem.
Solving the matrix eigenvalue equation amounts to matrix diagonalization. If you give your matrix H' to a diagonalization routine it will return a vector e with eigenvalues and a matrix V whith the corresponding eigenvectors as columns, such that \[ H'_{ij}V_{jk} = e_k V_{ik} \;. \]
gsl_matrix structure.
Hint:
int n=20;
double s=1.0/(n+1);
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/s/s);
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);
gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_VAL_ASC);
fprintf (stderr, "i exact calculated\n");
for (int k=0; k < n/3; k++){
double exact = pi*pi*(k+1)*(k+1);
double calculated = gsl_vector_get(eval,k);
fprintf (stderr, "%i %g %g\n", k, exact, calculated);
}
for(int k=0;k<3;k++){
printf("%g %g\n",0.0,0.0);
for(int i=0;i<n;i++) printf("%g %g\n",(i+1.0)/(n+1),gsl_matrix_get(evec,i,k));
printf("%g %g\n",1.0,0.0);
printf("\n\n");
}