program harmonic0 * * Risolve l'oscillatore armonico quantistico. * Integrazione solo in avanti con metodo di Numerov. * Ricerca degli autovalori mediante "shooting method". * implicit none * Numero massimo di punti nella griglia integer mshx parameter (mshx = 2000) real*8 pi parameter (pi=3.14159265358979323846d0) integer mesh, i, icl integer nodes, hnodes, ncross, kkk, iterate real*8 xmax, dx, ddx12, xmcl, norm real*8 eup, elw, e real*8 x(0:mshx), y(0:mshx), p(0:mshx), vpot(0:mshx), f(0:mshx) character*80 fileout * * Input e inizializzazione della griglia * print '($,a)','Massimo x (valore tipico: 10) ? ' read*, xmax print '($,a,i5,a)', $ 'Numero punti griglia (massimo:',mshx,') ? ' read*,mesh if (mesh.gt.mshx) stop dx = xmax/mesh ddx12 = dx*dx/12.d0 * * Definizione del potenziale (simmetrico rispetto a x=0) * do i = 0, mesh x(i) = dfloat(i) * dx vpot(i) = 0.5d0 * x(i)*x(i) enddo print '($,a)','Nome file output = ' read '(a)',fileout open(1,file=fileout,form="formatted",status="new") * * Entrata loop di ricerca nuovo autovalore * 999 continue * * Input numero nodi (<0 per terminare) * print '($,a)', 'Numero nodi (-1=esci) ? ' read*,nodes if (nodes.lt.0) then close(1) stop endif * * Definizione limiti all'autovalore * eup = vpot(mesh) elw = eup do i=0,mesh elw = min(elw,vpot(i)) eup = max(eup,vpot(i)) enddo * * Definizione energia per singola valutazione * print '($,a)','Energia da provare (0=ricerca con bisezione) ? ' read*,e if ( e .eq. 0.d0 ) then * Modalita' ricerca automatica autovalore e = 0.5d0 * (elw + eup) iterate = 1 else * Modalita' test di un singolo autovalore iterate = 0 endif kkk = 0 * * Punto di entrata per calcolo soluzione a energia e * 1 continue kkk = kkk + 1 * * Definizione funzione g usata nell'algoritmo di Numerov * f(0) = ddx12*(2.d0*(vpot(0)-e)) do i=1,mesh * f <= 0 regione permessa classicamente * f > 0 regione proibita classicamente f(i)=ddx12*2.d0*(vpot(i)-e) * Nessun f(i) puo' essere zero altrimenti * il cambio di segno non viene rilevato. * La linea seguente e' un trucco per evitare * di imbattersi in f(i) nulli: if( f(i) .eq. 0.d0 ) f(i)=1.d-20 * Registra l'indice a cui e' avvenuto * l'ultimo cambio di segno: if( f(i) .ne. sign(f(i),f(i-1)) ) icl=i enddo if(icl.ge.mesh-2) then print*,'Ultimo cambio segno troppo lontano.' stop endif * * f richiesta dal metodo di Numerov: * do i=0,mesh f(i) = 1.0d0 - f(i) enddo do i=0,mesh y(i) = 0.d0 enddo * * Determinazione funzione d'onda nei primi due punti * hnodes = nodes/2 if (2*hnodes.eq.nodes) then y(0) = 1.0d0 y(1) = 0.5d0*(12.d0-10.d0*f(0))*y(0)/f(1) else y(0) = 0.d0 y(1) = dx endif * * Integrazione verso l'esterno e conteggio dei * cambi di segno * ncross=0 do i=1,mesh-1 y(i+1) = ((12.d0-10.d0*f(i))*y(i)-f(i-1)*y(i-1))/f(i+1) if ( y(i) .ne. sign(y(i),y(i+1)) ) ncross=ncross+1 enddo * * Verifica il numero di cambi di segno * print '(2i4,f14.8)',kkk, ncross, e if ( iterate .ne. 0 ) then if (ncross.gt.hnodes) then * Avvenuti piu' cambi di segno del previsto. * L'energia e' troppo alta e va abbassata * (--> semiintervallo inferiore). eup = e else * Cambi di segno in numero giusto o inferiore. * L'energia e' troppo bassa e va alzata * (--> semiintervallo superiore). elw = e endif * Nuovo valore da provare: e = 0.5d0 * (eup+elw) * Criterio di convergenza: if (eup-elw .gt. 1.d-10) go to 1 endif * * La convergenza e' stata raggiunta, o non erano richieste * iterazioni. Attenzione! y e' una funzione d'onda non * normalizzata. * * Rinormalizzazione della funzione d'onda in modo che * int |psi|^2 dx = 1: *** NON LO FACCIAMO *** * Il problema e' la divergenza a grandi x. * * norm = 0.d0 * do i=1,mesh * norm = norm + y(i)*y(i) * enddo * norm = dx * (2.d0 * norm + y(0)*y(0)) * norm = sqrt(norm) * do i=0,mesh * y(i) = y(i) / norm * enddo * * Calcolo della densita' classica relativa all'energia e * (per confronto): * xmcl = sqrt(2.d0*e) do i=icl,mesh p(i) = 0.d0 enddo do i=0,icl-1 p(i) = 1.d0/sqrt(xmcl**2 - x(i)**2)/pi enddo * Semiasse negativo: do i=mesh,1,-1 write (1,'(f7.3,3d16.8,f12.6)') $ -x(i),(-1)**nodes*y(i),y(i)*y(i),p(i),vpot(i) enddo * Semiasse positivo: do i=0,mesh write (1,'(f7.3,3d16.8,f12.6)') % x(i),y(i),y(i)*y(i),p(i),vpot(i) enddo go to 999 end