Applicazione della teoria perturbativa

Particolare attenzione merita la sezione del programma:

 i = icl
 ycusp = (y(i-1)*f(i-1)+f(i+1)*y(i+1)+10.d0*f(i)*y(i)) / 12.d0
 dfcusp = f(i)*(y(i)/ycusp - 1.d0)
 *  Aggiornamento autovalore usando teoria delle perturbazioni
 de = dfcusp/ddx12 * ycusp*ycusp * dx

il cui scopo è quello di stimare al primo ordine in teoria delle perturbazioni la differenza $\delta e$ fra l'autovalore attuale e quello esatto.

Vogliamo spiegare il meccanismo che sta alla base di questa stima (prima di affrontare questa sezione occorre essere familiari col metodo perturbativo descritto nel prossimo capitolo). Ricordiamo che icl è l'indice corrispondente al punto di inversione classico. L'integrazione viene effettuata in avanti sino a questo indice, ed all'indietro pure sino a questo indice. icl è quindi l'indice di raccordo tra le due funzioni. La funzione di destra viene riscalata in modo che non vi sia una discontinuità della funzione $y$ al punto di raccordo, tuttavia la derivata prima $dy/dx$ sarà in generale discontinua (a meno che l'autovalore non sia quello cercato).

Si può quindi dire che y(icl) è il valore previsto dalla ricorrenza di Numerov calcolata usando come punto centrale sia icl-1 che icl+1. La ricorrenza di Numerov usando come punto centrale icl non è però mai stata applicata. Il valore calcolato per ycusp è appunto quanto predice la ricorrenza di Numerov usando icl come punto centrale; il problema è che è diverso da y(icl).

L'ottica in cui ci porremo ora è quella di pensare che la funzione ottenuta sia la soluzione esatta di un problema diverso, e precisamente di un problema in cui il potenziale nel solo intervallo $\Delta x$ centrato su y(icl) è diverso da quello dato. Il potenziale modificato ``curva'' la soluzione nel modo osservato. Una volta trovata la modifica del potenziale necessaria per ottenere questo effetto, la teoria delle perturbazioni ci fornisce una stima della differenza di autovalore rispetto all'originale.

Riesaminiamo la formula di Numerov (2.44), e notiamo come la formula (pensiamola applicata al caso n=icl-1) ci fornisce in realtà il solo prodotto y(icl)f(icl). Solitamente da questo prodotto estraiamo y(icl) perchè f(icl) si pensa dato. Supponiamo ora invece che f al punto icl abbia un diverso valore fcusp non noto, ma tale che la funzione trovata soddisfi alla relazione di Numerov anche nel punto icl. Deve allora essere

               fcusp ycusp = f(icl) y(icl)

in quanto questo prodotto è quanto fornito dal metodo di Numerov (integrando da icl-1 o icl+1), e ycusp è il valore che y deve avere affinchè la ricorrenza di Numerov sia soddisfatta anche in icl. Questo spiega la definizione di dfcusp (variazione di fcusp) calcolata dal programma.

Il prossimo step è quello di calcolare la corrispondente variazione del potenziale. Per ``potenziale'' intendiamo qui il potenziale unidimensionale efficace che appare nell'equazione (4.53), $W(x)=r^2 V + (\ell+1/2)^2$, e da qui nelle definizioni per gli f(i). Differenziando si trova che la variazione del potenziale $\delta W$ che dà luogo a una variazione $\delta f$ è data da

$$\delta W = - \frac{12}{(\Delta x)^2} \delta f$$ (4.54)

La teoria delle perturbazioni al primo ordine ci dà allora la corrispondente variazione dell'autovalore:

$$\delta e = \langle y\vert\delta W\vert y\rangle = \int \vert... ...ert^2 \delta W (x) dx = \vert y(x_c)\vert^2 \delta W \Delta x$$ (4.55)

avendo sviluppato l'integrale in una sommatoria sulla griglia, dove l'unico contributo non nullo viene fornito dalla regione di ampiezza $\Delta x$ centrata sul punto di inversione classico $x_c$ (ossia, nel programma, l'indice icl). Abbiamo dunque trovato che

$$\delta e = - \frac{12}{(\Delta x)^2} \vert y(x_c)\vert^2 \delta f \Delta x$$ (4.56)

è la differenza di autovalore per passare dal potenziale originale a quello perturbato. A parte il diverso segno--in quanto noi desideriamo risalire all'autovalore del problema non perturbato--questa è esattamente l'espressione usata dal programma per il calcolo della correzione all'autovalore de.

Finchè si è lontani dalla convergenza, questa correzione è grossolana e può portare a errori considerevoli. Per questo motivo il programma è ``protetto'' dalla linea

 e = max(min(e+de,eup),elw)

che impedisce in tutti i casi che venga stimato un nuovo valore di e all'esterno dell'intervallo [elw,eip]. Quando il programma procede verso la convergenza, la stima diventa via via sempre migliore e consente una convergenza assai rapida nella fase finale.