Resolución numérica de la ecuación de Schrodinger unidimensional.


Se resolvió el problema de una partícula de masa igual a la del electrón a velocidades no relativistas, utilizando lenguaje C y exportando los resultados en formato texto para poder graficar los resultados mediante Octave.

La idea de la publicación, es ser una base para proyectos de este tipo ya que es algo difícil de encontrar en internet. El hecho de que sea resuelto en Lenguaje C (compilado con GCC) hace este metodo aún mas raro. y tiene la ventaja de aprovechar a fondo el poder de calculo disponible (a diferencia de utilizar Octave o Matlab).

El código fuente se puede ver a continuación:



/*
Autor: Fernando Filippetti - #93684
extracto de codigo parte de un trabajo práctico para la materia FISICA 3 - FIUBA.
Código fuente bajo licencia GNU GPLv3.
 
stdio.h
stdlib.h
math.h 
time.h
*/ 
 
#define offset  (20e-9)
#define pi  (3.141592653589793)
#define m  (9.1e-31)
#define c  (1.6e-19)
#define h  (6.62606957e-34)
#define hb  (h/(2*pi))
#define delta_t (0.0001e-15)
#define delta_x (0.1e-9)      //Paso minimo (0.1nmetro)
#define tmax  (450e-15)     //450femtosegundpos
#define xmax  ((2*offset)+150e-9)
#define X_0  (offset+10e-9)        // xo = 10 nanometros
#define sigma  (3.51e-9)     //varianza = 3.5nm
#define v  (2.085e5)
#define lambda  (h/(m*v))
#define ra1  ((hb*delta_t)/(2*m*pow(delta_x,2)))
#define tamx ((uint)(xmax/delta_x))
#define tamt ((uint)(tmax/delta_t))
#define T0  0
#define T150 (int)(tamt*(151.0/450))
#define T250 (int)(tamt*(251.0/450))
#define T350 (int)(tamt*(351.0/450))
#define T450  tamt
#define POT_BARR (0.025*1.6e-19)


typedef double presicion_t;

int main()
{
FILE *fp1,*fp2,*fp3,*fp4,*fp5;
fp1 = fopen("T0.txt","w");
fp2 = fopen("T1.txt","w");
fp3 = fopen("T2.txt","w");
fp4 = fopen("T3.txt","w");
fp5 = fopen("T4.txt","w");
clock_t start = clock();

int i,t=0,it,ix=0;
fprintf(stderr,"\t Calculador de ecuacion de onda. \n\n\t t=%d x=%d ra=%f \n\t Tmax:%1.2ffs \t Xmax:%1.2fnm \n\n",tamt,tamx,ra1,tmax/1e-15,xmax/1e-9);
// double_t T[(int)(tmax/delta_t)]; ,j,k
// double_t X[(int)(xmax/delta_x)];
// double phir[tamx][tamt],phii[tamx][tamt];

presicion_t** phir;
phir = (presicion_t**) malloc(2*sizeof(presicion_t*));
phir[t] = (presicion_t *) malloc(tamx*sizeof(presicion_t)); // pido memoria para t, y para t+1 
phir[t+1] = (presicion_t *) malloc(tamx*sizeof(presicion_t));

presicion_t** phii;
phii = (presicion_t**) malloc(2*sizeof(presicion_t*));
phii[t] = (presicion_t *) malloc(tamx*sizeof(presicion_t));
phii[t+1] = (presicion_t *) malloc(tamx*sizeof(presicion_t));

presicion_t V[tamx];
presicion_t sum=0,ra; 
ra=ra1; // calculo "ra" una sola vez y despues lo uso (para no hacer POW() tantas veces.

/**************************************************************************************************
for(i=1;i!=tamx;i++){ //inicializo el potencial (cuadrado)
 if(offset+70e-9 < i*delta_x && i*delta_x <= offset+ 75e-9)
  {V[i]= -0.025*1.6e-19;} //energia del pozo de potencial en joules.
 else {V[i]=0;}
}*/
/**************************************************************************************************/
for(i=1;i!=tamx;i++){//inicializo el potencial (trianglar)
 if(offset+70e-9 < i*delta_x && i*delta_x <= offset+ 72.5e-9)
  {V[i]= ((i*delta_x-(offset+70e-9))/2.5e-9) * POT_BARR;//energia del pozo de potencial en joules.
  } 
 else if(offset+72.5e-9 < i*delta_x && i*delta_x <= offset+75e-9 )
  {V[i]= (1-((i*delta_x-(offset+72.5e-9))/2.5e-9)) * POT_BARR; 
  }
 else {V[i]=0;}
}
/**************************************************************************************************/
for (i=0; i!=tamx; i++) {//inicializo la primera fila de cada matriz con la funcion de onda inicial.
    phir[0][i]= cos(2*pi*(((i*delta_x)-X_0)/lambda)) *exp(-pow(((i*delta_x)-X_0)/sigma,2));
        phii[0][i]= sin(2*pi*(((i*delta_x)-X_0)/lambda)) *exp(-pow(((i*delta_x)-X_0)/sigma,2));
 sum+=pow(phir[0][i],2)+pow(phii[0][i],2);
    }
fprintf(stderr,"\t Coeficiente de normalizacion:%f \n",1/sqrt(sum));
/**************************************************************************************************/
sum=sqrt(sum);
for (i=0; i!=tamx; i++){ // ¡normalizo la cuacion de onda inicial!
         phir[0][i]= (phir[0][i])/sum;
         phii[0][i]= (phii[0][i])/sum;
     }
/**************************************************************************************************/
for (ix=0; ix != tamx; ix++){ //imprimo al archivo la ecuacion de onda inicial
  fprintf(fp1,"%e \t, %e \t, %e , %e , \n",(ix*delta_x)-offset,phir[T0][ix],phii[T0][ix],V[ix]/c);
  }
/**************************************************************************************************/
for (it=0; it != tamt+1; it++){ //resuelvo la ecuacion diferencial propiamente dicha, de manera iterativa (de otra forma harian falta mas de 32Gb de RAM.
  for (ix=1; ix != tamx; ix++){
         phir[(it+1)%2][ix]=phir[(it%2)][ix]-(ra*(phii[(it%2)][ix-1]-2*phii[(it%2)][ix]+phii[(it%2)][ix+1]))+(delta_t/hb)*V[ix]*phii[(it%2)][ix];
         phii[(it+1)%2][ix]=phii[(it%2)][ix]+(ra*(phir[(it%2)][ix-1]-2*phir[(it%2)][ix]+phir[(it%2)][ix+1]))-(delta_t/hb)*V[ix]*phir[(it%2)][ix];
             #ifdef NORM 
  sum+=phir[(it+1)%2][ix]*phir[(it+1)%2][ix]+phii[(it+1)%2][ix]*phii[(it+1)%2][ix]; 
  #endif
  }
  #ifdef NORM /* Esta funcion se habilita definiendo NORM*/
  sum=sqrt(sum);
  for (i=0; i !=tamx; i++){   //normalizo cada funcion de onda generada
         phir[(it+1)%2][i]= (phir[(it+1)%2][i])/(sum);
         phii[(it+1)%2][i]= (phii[(it+1)%2][i])/(sum);
      } 
  #endif
 if(it==T150){ for (ix=0; ix != tamx; ix++){
   fprintf(fp2," %e \t, %e \n",phir[(it+1)%2][ix],phii[(it+1)%2][ix]);
   }
  }
 else if(it==T250){for (ix=0; ix != tamx; ix++){
   fprintf(fp3," %e \t, %e \n",phir[(it+1)%2][ix],phii[(it+1)%2][ix]);
   }
  }
 else if(it==T350){for (ix=0; ix != tamx; ix++)        {
   fprintf(fp4," %e \t, %e \n",phir[(it+1)%2][ix],phii[(it+1)%2][ix]);
   }
  }
 else if(it==T450){for (ix=0; ix != tamx; ix++){
   fprintf(fp5," %e \t, %e \n",phir[(it+1)%2][ix],phii[(it+1)%2][ix]);
   }
  }
    }
/**************************************************************************************************/
fclose(fp1);// cierro los archivos y libero la memoria pedida al sistema.
fclose(fp2);fclose(fp3);fclose(fp4);fclose(fp5);
for(i=0; i<2; i++){
free(phir[i]);
free(phii[i]);}
free(phir); free(phii);
fprintf(stderr,"\t Tiempo de computo: %ld Minutos.   (%lds.)\n",((int)clock()-start)/(60*CLOCKS_PER_SEC),((int)clock()-start)/CLOCKS_PER_SEC);
return 0;
}  


Con lo cual despues de un par de minutos se pueden graficar los resultados que quedan de esta forma:





Comentarios

Entradas populares