Programma: harmonic1

harmonic1.c (versione C) o harmonic1.f (versione Fortran 77) costituisce una versione migliorata di harmonic0, che evita il problema della divergenza a grandi $x$.

Il trucco utilizzato per evitare la divergenza è quello di effettuare due integrazioni: una in avanti partendo da $x=0$, e una all'indietro partendo da $x_{\rm max}$. L'autovalore viene fissato dalla condizione che il raccordo fra le due parti sia continuo con derivata prima continua.

Il punto di raccordo scelto coincide col punto di inversione classica $x_{\rm cl}$, tale che $V(x_{\rm cl})=E$. Tale punto dunque varia man mano che diversi $E$ vengono provati, ed è ovviamente determinato entro la precisione che può fornire la griglia. In pratica viene identificato l'indice icl del primo punto della griglia $x_c = {\tt icl}\Delta x$ tale che $V(x_c)>E$. Possiamo quindi solo dire che il punto di inversione classico è situato fra $x_c-\Delta x$ e $x_c$.

L'integrazione verso destra viene portata avanti fino a icl, ottenendo così una funzione $\psi_L(x)$ definita in $[0,x_c]$, contando ancora il numero di cambiamenti di segno $n$ come in harmonic0. Se $n$ non è quello corretto l'energia è abbastanza distante dal valore richiesto, e si procede ad un aggiustamento per bisezione come descritto in precedenza (senza che sia necessario integrare al di là di $x_c$, regione in cui si sa a priori che non vi possono essere nodi). Se invece il numero di nodi è quello giusto, il programma procede con l'integrazione da destra verso sinistra^2.1, fermandosi allo stesso indice icl corrispondente a $x_c$ e ottenendo così una funzione $\psi_R(x)$ definita in $[x_c,x_{\rm max}]$.

A quel punto vi sono due valori della funzione d'onda in $x_c$: $\psi_L(x_c)$ e $\psi_R(x_c)$. La prima operazione effettuata è quella di riscalare $\psi_R(x)$ di un fattore $\psi_L(x_c)/\psi_R(x_c)$, in modo che via sia un raccordo continuo fra le due funzioni in $x_c$. Fatto questo, l'intera $\psi(x)$ viene rinormalizzata in modo che sia $\int \vert\psi(x)\vert^2 dx = 1$.

Ora inizia la parte nuova e cruciale: il calcolo della discontinuità della derivata prima alla giunzione, $\psi'_R(x_c) - \psi'_L(x_c)$. Questa differenza dovrebbe pure essere nulla per una buona soluzione, ma questo non sarà vero se non (all'interno di una precisione specificata) quando $E=E_n$. Il segno della differenza ci permette di capire se $E$ è troppo alta o troppo bassa, e quindi di applicare nuovamente il metodo di bisezione per migliorarne la stima.

Per calcolare la discontinuità, utilizziamo (indicando per brevità con $i$ l'indice corrispondente al punto di griglia $\tt icl$):

$$\begin{array}{l} y^L_{i-1} = y^L_i - y'^L_i \Delta x + \frac... ... \frac{1}{2} y''^R_i (\Delta x)^2 + O[(\Delta x)^3] \end{array}$$ (2.45)

Notiamo che $y^L_i=y^R_i=y_i$, e anche $y''^L_i=y''^R_i=y''_i=-g_i y_i$, come garantisce il metodo di Numerov. Sommando fra loro le due equazioni abbiamo allora

$$y^L_{i-1} + y^R_{i+1} = 2y_i + ( y'^R_i - y'^L_i ) \Delta x -g_i y_i (\Delta x)^2 + O[(\Delta x)^3]$$ (2.46)

da cui

$$y'^R_i - y'^L_i = \frac { y^L_{i-1} + y^R_{i+1} - [ 2 - g_i (\Delta x)^2 ] y_i } { \Delta x } + O[(\Delta x)^2]$$ (2.47)

o anche, utilizzando la notazione (2.42),

$$y'^R_i - y'^L_i = \frac { y^L_{i-1} + y^R_{i+1} - ( 14 - 12f_i ) y_i } { \Delta x } + O[(\Delta x)^2]$$ (2.48)

Il programma calcola in questo modo la discontinuità nella derivata prima. Se il suo segno è positivo, l'energia è troppo alta e viene quindi scelto il semiintervallo inferiore; se negativo, l'energia è troppo bassa e viene scelto il semiintervallo superiore. Al solito, si ritiene di essere giunti a convergenza quando l'ampiezza dell'intervallo di energie è diventata minore di una certa piccola tolleranza.

Oltre alla funzione d'onda finale scritta sul file di output in modo analogo ad harmonic0 (ma ora è normalizzata!), il programma emette una linea sullo standard output per ogni iterazione, contenente quattro numeri: il numero dell'iterazione, il numero di nodi trovati (sul solo semiasse degli $x$ positivi), la discontinuità nella derivata $y'^R_i - y'^L_i$ (solo quando è stata calcolata, ossia se il numero di nodi era quello giusto--negli altri casi viene riportato zero) e la stima corrente dell'autovalore dell'energia.