\documentstyle[12pt]{article} \textwidth 16cm \textheight 23cm \topmargin 0cm \evensidemargin 0cm \oddsidemargin 0cm \parindent=0pt \parskip=10pt \begin{document} \begin{center} \Large \bf Numerieke methoden in de natuurwetenschappen\\ \large \end{center} \vspace{0.4cm} \noindent {\bf Ex.~3a: Bound states in a finite square well, solution by use of NAG library}\\ Consider the problem of calculating the bound states in a potential well of finite height you have already solved. Solve it again by calling one of the routines of the NAG library for the root searching. The NAG library is a collection of ready-to-use subroutines which can be called from a f77, f90 or C program. It is useful to learn using them. The first step is to identify the appropriate subroutine, the second to understand which is the required input and where the output is given. A list of the avalaible subroutines can be found by invoking nagdtext or by looking at the manuals on paper which can be found in the computer room of mathematics on the second floor. A description of how to use them can be found at: http://www.sci.kun.nl/cncz/software/local/html/nag.html. In short, you have to give the path by typing: setenv LD\_LIBRARY\_PATH /usr/lib:/usr/ccs/lib:/usr/local/lib and compile with the options: f90 myprog.f90~~ -lnag -lF77 or (hopefully, let me know if it works) CC myprog.c %-L/opt/SUNWspro/lib -R/opt/SUNWspro/lib -lnag -lF77 -lM77 -lsunmath -lm -lc The -lF77 option is needed because the NAG library is written in F77. Find a subroutine which performs a root searching and insert a call to it in your program. Check the results against your own. \end{document} Consider a square potential well of finite height V$_0$ and width $2a$, as in figure ~\ref{fig:put}. \begin{figure}[hbt] \centering \setlength{\unitlength}{0.00050000in}% % \begin{picture}(6066,5493)(2968,-6700) \thicklines \put(6001,-6361){\line( 0, 1){4800}} \put(3001,-2761){\line( 1, 0){1200}} \put(4201,-2761){\line( 0,-1){3600}} \put(4201,-6361){\line( 1, 0){3600}} \put(7801,-6361){\line( 0, 1){3600}} \put(7801,-2761){\line( 1, 0){1200}} \put(6001,-1561){\line( 0, 1){ 75}} \put(6001,-1486){\vector( 0, 1){ 75}} \put(6001,-6661){\small 0} \put(3901,-6661){\large $-a$} \put(7801,-6661){\large $a$} \put(6226,-1411){\large $V(x)$} \put(6301,-2761){\large $V_{0}$} \end{picture} \caption[1]{Finite potential well} \label{fig:put} \end{figure} The Schr\"odinger equation: \[ \frac{d^2\Psi}{dx^2}-\frac{2m}{\hbar^2}\left[V(x)-E\right]\Psi=0 \] with $V(x)$ given by: \begin{eqnarray*} V(x)&=&V_0 \qquad x<-a \mbox{~~~~~~~~~~~~~~~~~~~~~region $I$~~~}\\ V(x)&=&0 \qquad -aa \mbox{~~~~~~~~~~~~~~~~~~~~~~~~region $III$~} \end{eqnarray*} has general solutions for the bound states at energy less than $V_0$ in the three regions given by: \begin{eqnarray*} \Psi_{I}(x)&=&A~\exp^{\beta x}+B~\exp^{-\beta x} \qquad x<-a\\ \Psi_{II}(x)&=&C~\sin\alpha x+D~\cos\alpha x \qquad -aa\\[1ex] \alpha&=&\sqrt{\frac{2mE}{\hbar^2}}\\ \beta&=&\sqrt{2m(V_0-E)/\hbar^2} \end{eqnarray*} The boundary conditions require continuity of the wavefunction and of its derivative at the boundary, whence \\[1ex] Even states: $B=E=0, C=0, D\neq0, A=F \qquad \alpha \tan \alpha a = \beta$\\ Odd states: $B=E=0 , C\neq0, D=0, A=-F \qquad \alpha \cot \alpha a = -\beta$\\ Therefore eigenvalues of the Schr\"odinger equation are found by solving numerically: \begin{center} \begin{tabular}{ccc} $\alpha tg(\alpha a)-\beta$=0 &for even states\\ $\alpha ctg(\alpha a)+\beta$=0 &for odd states \end{tabular} \end{center} A more convenient form is: \begin{center} \begin{tabular}{ccc} $\alpha sin(\alpha a)-\beta cos(\alpha a)$=0 &for even states\\ $\alpha cos(\alpha a)+\beta sin(\alpha a)$=0 &for odd states \end{tabular} \end{center} Use $2a$=10~\AA, $m=1~m_e$. It is useful to use $\hbar^2=7.6199682 m_e$~eV~\AA$^2$ \begin{enumerate} \item The repetitive nature of the tangent suggests that there might be several solutions, one per each cycle of the function. Establish the intervals in which roots exist. Remember that, in one dimension, a square well has always at least one bound state at $E\neq 0$. \item Find all bound states with precision 10$^{-4}$~eV by use of the bisection, Newton and secant methods for several values of $V_0$ between $V_0$=0.2~eV and $V_0$=5~eV. \item Make a plot of the eigenvalues as a function of $V_0$, comparing also to the limiting case of an infinite well: \begin{equation} E_n= {\hbar^2\over 2m_e}({n\pi\over 2a})^2. \end{equation} Do not forget to label the axis giving the used units and for which values of the parameters (e.g which value of $a$) the calculation is done. The idea is that, even without the text of the exercise, by reading the report one can understand and reproduce the calculation. Comment briefly the results. If you use gnuplot, you could plot the eigenvalues with linespoints so that you can trace back for which values of $V_0$ the calculations has been performed. \item Facultative: Make a plot of the eigenvalues as a function of $a$ for $V_0$=3 eV. %calculate the wavefunction corresponding to the first bound %state for two values of $V_0$ and comment the results. %\item Could you find out which is the value of $V_0$ which gives the first 3 %bound states at energy $E_1$=0.26230 $E_2$=1.03633 $E_3$= 2.27006 eV? \end{enumerate} \footnote{ More details on this problem and on computational strategy can be found in: P.L. De Vries {\it A first course in computational Physics}, J. Wiley \& sons, N.Y. 1994} \end{document}