Projet informatique C++ : Résolution équation de Schrodinger

Version imprimable

Voir 40 message(s) de cette discussion en une page

02/12/2014, 20h00
Andromedae

Merci beaucoup pour ta réponse !

Je vais continuer à me plonger dedans avec la notice explicative en parallèle du code de mon prof.
Je le vois demain donc je vais lui en toucher 2-3 mots histoire de mieux cerner le projet.

Bon j'ai pu modifier mon code, et il fonctionne plutot bien pour le moment.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include <iostream>
#include <cmath>
 
using namespace std;
 
#include <cstdlib>
#include <vector>
/*#include "matprod.h"*/
#include <fstream>
 
#include <stdexcept>
 
 
 
 
// -------------------------------------------------- //
// ------ Définition du potentiel carré simple ------ //
// -------------------------------------------------- //
 
const double PLANCK_CONSTANT = 1;
const double ELECTRON_MASSE = 1;
const double E = 1;
 
double V(double x, double xmin, double xmax, double V0)
{
	if (xmin < x && x < xmax) {
 
		return V0;
	}
	throw std::out_of_range(" x est defini hors du puit de potentiel ! Veuillez saisir x entre xmin et xmax ");
}
 
double SimpleSquarePotential(double E, double V0)
{
	double J = ((2 * ELECTRON_MASSE) / (PLANCK_CONSTANT*PLANCK_CONSTANT))*(E - V0);
 
	return J;
}
 
 
// -------------------------------------------------- //
// --------- Partie principale du programme --------- //
// -------------------------------------------------- //
 
int main()
{
	try {
 
		double v = V(0, 0, 2, -2);
		double j = SimpleSquarePotential(E, v);
		cout << j << endl;
 
	}
	catch (std::exception const & e) {
		std::cerr << e.what();
		return -1;
	}
 
 
	return 0;
}

Maintenant il ne reste plus qu'a s'attaquer au gros du projet ...

Puisqu'on parle de système d'unité, le truc de pro c'est de le faire façon CLHEP (Soit en incluant directement la lib, soit en ré-implémentant les deux fichiers qui vont bien)

http://proj-clhep.web.cern.ch/proj-c...its/units.html

On a alors
Code:

1 2 static const HepDouble h_Planck = 6.6260755e-34 * joule*s;
avec
Code:

1 2 3 4 5 static const HepDouble megaelectronvolt = 1. ; static const HepDouble electronvolt = 1.e-6*megaelectronvolt; static const HepDouble e_SI = 1.60217733e-19; // positron charge in coulomb static const HepDouble joule = electronvolt/e_SI;
Il suffit alors de diviser le résultat que tu obtiens par les unité dans lesquelles tu veux le présenter.

Mais c'est surement overkill pour un petit projet comme le tiens

Merci pour cette précision ;)
Je vais prendre l'unité atomique dans mon cas qui me simplifieras bien mon problème.

D'ailleurs je viens de rajouter mes 2 autres types de potentiels :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include <iostream>
#include <cmath>
 
using namespace std;
 
#include <cstdlib>
#include <vector>
#include "matprod.h"
#include <fstream>
 
#include <stdexcept>
 
 
// -------------------------------------------------- //
// ------ Définition des valeurs & constantes ------- //
// -------------------------------------------------- //
 
const double REDUCED_PLANCK_CONSTANT = 1;		// Définition de la constante de Planck réduite (unité atomique)
const double ATOMIC_MASSE = 1;				// Définition de la masse atomique
const double E = -1;					// Définition de l'énergie atomique en hartree (unité atomique)
double x ;
const double xmin = 0 ;
const double xmax = 2 ;
double Delta = (xmax-xmin)/2 ;
double Sigma = (xmax+xmin)/2 ;
double V ;
 
 
// -------------------------------------------------- //
// ------ Définition du potentiel carré simple ------ //
// -------------------------------------------------- //
 
double SimpleSquarePotential (double x, double xmin, double xmax, double V0)
{
	if (xmin < x && x < xmax) {
 
        return V = V0 ;
    }
    throw std::out_of_range(" V0 vaut l'infini ");
}
 
// -------------------------------------------------- //
// ------ Définition du potentiel quadratique ------- //
// -------------------------------------------------- //
 
double QuadraticPotential (double x, double xmin, double xmax, double V0, double Sigma, double Delta)
{ 
	if (xmin < x && x < xmax ) {
 
	double V = V0*((1-((x-Sigma)*(x-Sigma)))/(Delta*Delta)) ;
 
	return V ; }
 
	else 
 
	return V = 0 ;
}
 
// ------------------------------------------------------- //
// ------ Définition du potentiel de Lennard-Jones ------- //
// ------------------------------------------------------- //
 
double LennardJonesPotential (double x, double xmin, double xmax, double V0)
{
	double V = V0 * ( (pow((xmin/x),6)) - (pow((xmin/x),12)) ) ;
 
	return V ;
}
 
// --------------------------------------------------------------------- //
// ------ Définition du terme J(x) dans l'équation de Schrödinger ------ //
// --------------------------------------------------------------------- //
 
double J (double E, double V)
{
	double J = ((2 * ATOMIC_MASSE) / (REDUCED_PLANCK_CONSTANT*REDUCED_PLANCK_CONSTANT))*(E - V);
 
    return J;
}
 
// -------------------------------------------------- //
// --------- Partie principale du programme --------- //
// -------------------------------------------------- //
 
int main()
{
 
cout <<" ------------------------------------------------------------------ "<< endl;
cout <<" ------------------------------------------------------------------ "<< endl;
cout <<" ---------------------- Projet Informatique ----------------------- "<< endl;
cout <<" --- Résolution de l'équation de Schrödinger stationnaire à 1 D --- "<< endl;
cout <<" ------------------------------------------------------------------ "<< endl;
cout <<" ------------------------------------------------------------------ "<< endl;
 
 
// ---------------------------------------------- //
// ------ Calcul du potentiel carré simple ------ //
// ---------------------------------------------- //
 
    try {
 
        double SSP = SimpleSquarePotential(1, xmin, xmax, -18);
        double j1 = J(E, SSP);
        cout << "La valeur de J(x) pour le potentiel carré vaut : J(x) = " << j1 << endl ;
 
    }
    catch (std::exception const & e) {
        std::cerr << e.what();
        return -1;
    }
 
// ---------------------------------------------- //
// ------ Calcul du potentiel quadratique ------- //
// ---------------------------------------------- //
 
	double QP = QuadraticPotential (0,xmin,xmax,-18,Sigma,Delta) ;
	double j2 = J(E, QP) ;
		cout << "La valeur de J(x) pour le potentiel quadratique vaut : J(x) = " << j2 << endl;
 
// -------------------------------------------------- //
// ------ Calcul du potentiel de Lennard-Jones ------ //
// -------------------------------------------------- //
 
	double LJP = LennardJonesPotential (1,xmin,xmax,-18) ;
	double j3 = J(E, LJP) ;
		cout << "La valeur de J(x) pour le potentiel de Lennard-Jones vaut : J(x) = " << j3 << endl;
 
 
    return 0;
}

03/12/2014, 11h43
Andromedae

Maintenant l'objectif est de voir comment je peux définir mon équation de Schrodinger à l'ordre 2 en 2 équations à l'ordre 1.
03/12/2014, 12h50
Flob90

Tu poses une fonction comme étant la dérivée de ta fonction d'onde, et voila tu as ta décomposition en deux équations différentielles du premier ordre.
03/12/2014, 14h17
white_tentacle

Citation:

Envoyé par _zzyx_

Puisqu'on parle de système d'unité, le truc de pro c'est de le faire façon CLHEP (Soit en incluant directement la lib, soit en ré-implémentant les deux fichiers qui vont bien)

http://proj-clhep.web.cern.ch/proj-c...its/units.html

Ça me semble très loin de ce que permet un boost::unit. *Notamment, je n’ai rien vu qui empêche de mélanger allègrement des quantités de dimension différentes (et je passe sur les problème de précision / explosion numérique) que peuvent engendrer le fait que tout soit converti dans l’unité de base).

Bref, pas si pro que ça, et surtout, plutôt old-fashioned. On fait nettement mieux depuis.
09/12/2014, 12h35
Flob90

Je n'ai pas suffisamment compris comment traiter les conditions initiales pour la méthode. Comme je te l'ai dit dans un message précédent, personnellement j'irais demander à ton prof d'expliquer cette partie. De manière concrète, pour le premier potentielle que tu dois traiter, je suppose que les singularité sont au niveau des x_min et x_max, ça veut dire qu'il faut appliquer la méthode entre x_min+eps et X_max-eps ou une fois à x_min+eps et x_min+x_max/2 et une fois à x_min+x_max/2 et x_max-eps. Ensuite si on se place dans le premier cas, l'énoncé semble fixer (avec u la fonction d'onder) u(x_min)=u(x_max)=0 mais comment on choisit la valeur de départ de u'(x_min) et u'(x_max) pour appliquer la méthode de l'arbalète ?

Tout ça ce sont des questions qui n'ont rien à voir avec la programmation et ton prof sera surement une meilleur source de réponse que les forum.

Voir 40 message(s) de cette discussion en une page