Equation de schrodinger : calcul des états propres dans un puits (méthode type euler)

Soyez le premier à donner votre avis sur cette source.

Vue 11 436 fois - Téléchargée 425 fois

Description

Ce code code permet de calculer les etats propres d'une particule quantique (quanton) dans un puit de potentiel (symetrique et "infini").

Il se decompose en 4 etapes principales.

_Premiere etape: il vous demande d'entrer n (numero du mode que vous desirez) au clavier
_Deuxieme etape: il definit la parite de la fonction d'onde correspondante
_Troisieme etape: apres avoir fait diverses operations mineures il calcule le fondamental et tous les modes inferieurs au n que vous avez entrer (mais il ne les sauvegarde pas!), et affiiche l'energie correspondant a chaque mode.
_Quatrieme etape: il inscrit ensuite les valeurs de la fonction d'onde que vous avez demander dans un ficher nomme Energie.plt

Voila il ne vous reste plus qu'a admirer le resultal et tracant la fonction avec un logiciel comme gnuplot par exemple.

Le potentiel par defaut est un potentiel en x². Vous pouvez essayer diverses formes (puit carre ou bizaroides) ils marchent tous du moment qu'ils sont suffisament abrupts, profonds et symetriques (pour le changer faites le dans le source...).

Pour plus de detail sur son fonctionnement lisez les commantaires, j'ai essayer d'etre le plus clair possible.

Et si vous avez un doute faites un peu de theorie, ca n'a jamais tuer personne... Pour ceux que ca interressent les livres de LE BELLAC ("Physique quantique") et LEVY-LEBLOND ("Quantique") sont assez exaustifs pour une premiere approche (mais allez y doucement quand meme ;) )

Nb: La mise en page est assez moche si vous essayer de lire le texte avec un editeur classique, je vous conseille donc de le regarder sous vi (pour les puristes) ou xemacs (ou tout autre logiciel adapte a la programmation)

Source / Exemple :


/********************************************Inclusion des bibliotheques standard*******************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*********************************************Definition des constantes generales******************************************/

#define C 1.0 //Constante de Planck  Heu... ca serait cool de le prendre = a 1.... precision machine etc...
#define PI 3.141592654 //PI
#define l2 10 //Definition de la largeur du puit
#define h  0.001   //pas d'integration

/***********************************************Fonction definissant le potentiel*********************************************/

/*REMARQUE: plus le puit et profond plus je pourai resoudre de modes mais plus le programme sera lent.
 Exemple une profondeur de -100 permet de resoudre "seulement" les 12 premier modes du puit carre*/

double V (double i)
{
    return (-1000+(h*h*i*i/10));
}

main()
{

  /******************************************Toutes les variables dont je me sert******************************************/

 double *x;	                                         //Definition de la fonction d'onde
 double *v;	                                         //Definition de la derivee premiere de la fonction d'onde
 double r;                                           //Definition du pas dans la recherche de l'energie
 double E;                                           //Definition de l'energie
 double U;                                           //Autre definition de l'energie utilisee pour la dichotomie
 double e;                                           //Autre definition de l'energie pour la recherche grossiere
 double l; 	                                         //Definition de la longueur de l'espace d'integration
 int p,q,w;                                          //Sert lors de la recherche des modes
 double m; 	                                         //Masse de la particule
 int parite; 	                                     //Definit la parite de la fonction
 FILE *pf1,*pf2,*pf3;                                //Fichier dans lequel je vais stocker mes donnees
 int i,j,k,n,b,lg,z; 	                             //Definition de(s) truc(s) a incrementer dans les boucles
 

  /*******************************************Initialisation des variables**************************************************/

 printf("Quel mode cherchez vous? (0 correspond au fondamental)\nn=");
 scanf("%d",&n);                                            //n definit quel mode l'on recherche
 printf("\n\n");

 b=n;                                                       //On sauvegarde n qqpart car on va l'utiliser et le changer apres

 if (n%2==0)
    parite=1;                                               //On impose la parite de la fonction d'onde suivant le mode que l'on recherche
 else
     parite=0;
     

 

w=1;
 for (i=0;i<n;i++)
    {                                                        //On calcule le nombre de fois que la fonction va changer
      if (i%2==0)                                            //de facon de diverger
	w=w+2;
      else
	w=w+1; 
    }                                             

 

 l=20;     	                                                 //N sert a definir la longueur sur laquelle on integre
 lg=(int)(l/h)+1;                                            //Eje decoupe l'axe sur lequel j'integre en petits segments de longueur h
 m=1;      		                                             //Fmasse de la particule
 E=((C*C)/(2*m))*((PI/(2*l2))*(PI/(2*l2)))*(n+1)*(n+1)-1000; //RDefinition de l'energie pour le puits carre "infini"

                                                             //ARemarque si on chage la profondeur du puit de potentiel il faut changer e aussi

 e=-1000;                                                    //AEnergie minimale du puits
 r=0.01;                                                     //LPas dans la recherche de l'energie
 
 
 x=calloc(lg,sizeof(double));                                //Initialisation des tableaux
 v=calloc(lg,sizeof(double));
 

 
 //si la fonction est paire on lui impose sa valeur en zero comme etant 1
 //si ma fonction d'onde est paire la derivee est nulle en zero 
 //si la fonction est impaire on lui impose sa valeur en zero comme etant 0
 //si la fonction est impaire sa derivee est 1    
 if (parite == 1) 
   {
     x[0]=1.0;
     v[0]=0.0;     
   }
 else 
   {
     x[0]=0.0;        
     v[0]=1.0;        
   }
 
 
 
 pf1=fopen("Energie.plt","w");                               //Fichiers dans lesquels je vais preserver mes resultats
 
 
 n=0;                                                        //on initialise ces valeurs a zero a fin de faciliter
 p=0;                                                        //l automatisation du programme
 q=0;

 
 for(i=0;i<w;i++)
   {
     
     
     
     /****************************Recherche primaire de l'energie du niveau fondamental********************************/

     
     if (n%2==0)
       parite=1;                                              //On impose la parite de la fonction d'onde suivant le mode que l'on recherche
     else
       parite=0;
     

     
     if (parite == 1) 
       {
	 x[0]=1.0;
	 v[0]=0.0;                                               //Et l'on donne les conditions initialles correspondantes
       }
     else 
       {
	 x[0]=0.0;        
	 v[0]=1.0;        
       }
     

     if(i%2==0)
       {
	 for(e;e<-98;e+=r)
	   {
	     for(z=0;z<lg-1;z++)
	       {
		 v[z+1]=v[z]+h*(-((2*m)/(C*C)))*(e-V((double)z))*x[z];    //Algorithme d'integration: sorte de runge-kutta
		 x[z+1]=x[z]+h*v[z+1];                                    //a l'ordre 1
	       }
	     if   (x[(int)((l/h))] <0)                                //Si la facon de diverger change c'est que j'ai
	       {                                                      //depasser l'energie exacte du fondamental
		 E=e;                                                     //Ici diverger est traduit par:
		 U=E-r;                                                   //"si la valeur de ma fonction au dela du potentiel
		 break;                                                   //vaut moins qu'une certaine valeur"
	       }                                                      //U sert a preserver une valeur pour laquelle
	   }                                                          //la divergence est differente afin d'obtenir
	                                                              //les deux valeurs entre lesquelles je vais
	                                                              //faire ma dichotomie
	 
	 /************************Recherche dichotomique de l'energie du niveau fondamental********************************/
	 
	 if (n%2==0)
	   parite=1;                                                       //On impose la parite de la fonction d'onde suivant le mode que l'on recherche
	 else
	   parite=0;
	 
	 
	 if (parite == 1) 
	   {
	     x[0]=1.0;
	     v[0]=0.0;                                                     //Et l'on donne les conditions initialles correspondantes 
	   }
	 else 
	   {
	     x[0]=0.0;        
	     v[0]=1.0;        
	   }
	 
	 
	 for(j=0;j<100;j++)
	   {
	     for(k=0;k<lg-1;k++)
	       {
		 v[k+1]=v[k]+h*(-((2*m)/(C*C)))*(((E+U)/2)-V((double)k))*x[k];         //Je test la convergence pour une valeur 
		 x[k+1]=x[k]+h*v[k+1];                                                 //comprise entre E et U
	       }
	     
	     if   (x[(int)((l/h))] <0)                                             //Si cette convergence ce fait de maniere identique a celle de E trouvée
	                                                                           //precedement alors c'est une energie plus proche de la valeur exacte
	       E=(E+U)/2;
	     
	     if (x[(int)((l/h))]>0)                                                //En revenche si la convergence ce fait de manire differente alors
	                                                                           //je peux augmenter ma borne inferieure (U) car mon energie n'est
	       U=(E+U)/2;                                                          //pas dans le demi intervalle inferieur de [E;U]
	     
	     
	   }
       }
     /*REMARQUE: la divergence tantot vers plus l'infini tantot vers moins l'infini explique la presence de la
       deuxieme recherche primaire et dichotomique*/

     if(i%2 != 0)
       {

	 for(e;e<-98;e+=r)
	   {
	     for(z=0;z<lg-1;z++)
	       {
		 v[z+1]=v[z]+h*(-((2*m)/(C*C)))*(e-V((double)z))*x[z];        //Algorithme d'integration: sorte de runge-kutta
		 x[z+1]=x[z]+h*v[z+1];                                        //a l'ordre 1
	       }
	     if   (x[(int)((l/h))] >0)                                    //Si la facon de diverger change c'est que j'ai
	       {                                                          //depasser l'energie exacte du fondamental
		 E=e;                                                         //Ici diverger est traduit par:
		 U=E-r;	                                                      //"si la valeur de ma fonction au dela du potentiel
		 break;                                                       //vaut moins qu'une certaine valeur".
	       }                                                          //U sert a preserver une valeur pour laquelle
	   }                                                              //la divergence est differente afin d'obtenir
	                                                                  //les deux valeurs entre lesquelles je vais
	                                                                  //faire ma dichotomie
	 
	 
	 /************************Recherche dichotomique de l'energie du niveau fondamental********************************/
	 
	 for(j=0;j<100;j++)
	   {
	     for(k=0;k<lg-1;k++)
	       {
		 v[k+1]=v[k]+h*(-((2*m)/(C*C)))*(((E+U)/2)-V((double)k))*x[k];     //Je test la convergence pour une valeur 
		 x[k+1]=x[k]+h*v[k+1];                                             //comprise entre E et U
	       }
	     
	     if   (x[(int)((l/h))] >0)                                             //Si cette convergence ce fait de maniere identique a celle de E trouvée
	                                                                           //precedement alors c'est une energie plus proche de la valeur exacte
	       E=(E+U)/2;
	     
	     if (x[(int)((l/h))]<0)                                                //En revenche si la convergence ce fait de manire differente alors
	                                                                           //je peux augmenter ma borne inferieure (U) car mon energie n'est
	       U=(E+U)/2;                                                          //pas dans le demi intervalle inferieur de [E;U]	     
	   }

       }
     
     
     
     if (i==p )
       {
	 if (q%2==0)
	   {
	     p=p+2;
	     printf("Energie=%lf\n",E);                                            //Cette partie du code sert a distinguer les moment ou l'on doit
	   }                                                                       //rechercher une fonction paire de ceux ou l'on recherche une 
	 else                                                                      //fontion impaire (depend de la valeur de n).
	   {
	     p=p+1;
	     printf("Energie=%lf\n",E);
	   }
	 q++;
	 n=n+1;
       }
     
     

   }
 
 
 
 /**************************************Dessin de fonction d'onde correspondante*************************************/

 n=b;
 if (n%2==0)
   parite=1;                                                       //On impose la parite de la fonction d'onde suivant le mode que l'on recherche
 else
   parite=0;
 
 
 if (parite == 1) 
   {
     x[0]=1.0;
     v[0]=0.0;     
   }
 else 
   {
     x[0]=0.0;        
     v[0]=1.0;        
   }
 for(i=0;i<lg-1;i++)
   {
     fprintf(pf1,"%lf \t %lf \t %lf \n",h*i,x[i],V(i));            //On imprime les resultats dans un fichier et on complete suivant la paritee
     
     if (parite == 1)
       fprintf(pf1,"%lf \t %lf \t %lf \n",-h*i,x[i],V(i));
     
     if (parite != 1)
       fprintf(pf1,"%lf \t %lf \t %lf \n",-h*i,-x[i],V(i));
     
     
     v[i+1]=v[i]+h*(-((2*m)/(C*C)))*(E-V((double)i))*x[i];
     x[i+1]=x[i]+h*v[i+1];
     
     
     if (x[i]< -5)                                                //il ne sert a rien d'imprimer plus etant donne que une fois que la fonction a
       break;                                                     //commancer a diverger elle ne va pas s'arreter
     
     if (x[i] >5)
       break;
     
     
   }
 printf("\n\n");
 close(pf1);
}

Conclusion :


Je tiens a remercier specialement F.Hebert qui a coodonner ce projet.

J'essairai de donner une version plus complete de la facon dont marche ce programme ainsi que les ameliorations a y apporter ulterieurement.

Desole pour les maniaques de l'orthographe mais moi c'est pas mon truc alors....

Codes Sources

A voir également

Ajouter un commentaire Commentaires
Messages postés
1
Date d'inscription
jeudi 11 avril 2013
Statut
Membre
Dernière intervention
12 avril 2013

mercie nefraal c'est très interessant
Messages postés
3006
Date d'inscription
dimanche 14 avril 2002
Statut
Membre
Dernière intervention
31 décembre 2008

C'est mignon de croire au père noël à 22 ans.
Messages postés
1
Date d'inscription
vendredi 16 mars 2007
Statut
Membre
Dernière intervention
22 avril 2007

cher(e) ami(e),c'est exéllent ce que vous avez fait,je vous félicite,je désire que tu m'aide dans la resolution d'un puit fini et symetrique de largeure a et de profondeur vo>0,j'ai fait un premier programme mais ce qui me géne ce sont les unité de mesure donné à la machine ,j'èspere que tu m'envoie le programme quand tu finira a mon email: benalilou.h@hotmail.com pour que je puisse le comparer avec le mien est pourquoi pas se discuter dans les erreurs que j'ai mis de dont,merci au paravant;..
Messages postés
101
Date d'inscription
dimanche 1 mai 2005
Statut
Membre
Dernière intervention
13 décembre 2007
1
Est ce que tu pourrais mettre les fichiers source dans le zip?
Messages postés
3006
Date d'inscription
dimanche 14 avril 2002
Statut
Membre
Dernière intervention
31 décembre 2008

euhm, la partie code ça me bloque pas du tout, comme tu dis, c'est élémentaire. pas de souci de ce côté, non, c'est vraiment la méca quantique... j'ai aucune notion quantitative. j'ai lu qq bribes d'articles qualitatifs, mais...
Afficher les 8 commentaires

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.