Bonjour à tous.
Si je viens vous demander votre aide c'est que là je suis vraiment bloqué... Ca fait une bonne dizaine d'heures que je travail sur ce projet et ça avançait plutôt bien jusqu'à présent mais là je me retrouve dans une impasse.
Que je vous explique mon problème: Dans mon UE de méthodes numériques il m'est demandé de réaliser la modélisation de la fonction d'onde d'une particule à travers un potentiel (énoncé en annexe.).
La partie interessante de l'énoncé est la suivante:
Pour résoudre le problème numériquement, il est moins couteux en temps de calcul de prendre le problème à
l’envers :
– on prend en x > L une fonction d’onde eikx ;
– on écrit les conditions aux limites en L ;
– on résout l’équation différentielle entre L et 0 ;
– les conditions aux limites en 0 permettent de calculer la solution en x < 0 : Aeikx + Be−ikx ;
– enfin,A=1/T etB=R/T.
L'équation de Schrodinger dans un état stationnaire à une dimension est donné par l'équation différentielle de type y"=-a(x)y.
J'ai donc créé un code me permettent de résoudre cette équation avec a(x) variant suivant si on se trouve en dedans ou en dehors du potentiel.
Et....ca marche! (enfin ça m'a l'air). Petit graphe sous root en annexe.
Mais là où je bloque c'est pour déterminer le coefficient de transmission et le coefficient de réflexion. Suivant les conditions limite en x=0 on a psi(0)=A+B et psi(0)'=k(A-B). Donc on en tire A et B et donc R et T. Mais là je trouve des trucs complètement aberrants! Sachant que R+T=1 normalement. (mon code en annexe)
Est-ce que le ciel (ou vous) peut me venir en aide?! Merci.
PS: n'hésitez pas si vous avez des questions, j'ai conscience de ne pas avoir été forcément très clair sur tout les points^^
Cordialement, Jeremy.
Annexe:
Enoncé:
Graphe:
Code:
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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 #ifndef __EDO2_H__ #define __EDO2_H__ #include "Vecteur.h" #include <cmath> class EDO2 { protected: int N; double E; Vecteur yf; // Valeurs de y(xf) stockees ici par rk4() double *xtab; Vecteur *ytab, *yptab; void stocke(int i,double x, Vecteur y, Vecteur yp); public: double L=21e-9; double V=1.6e-19,m=1.9e-31,hb=1.04082e-33, k1=sqrt(2*m*E)/hb, k2=sqrt(2*m*(V-E))/hb; EDO2(int Ni, long double Ei); ~EDO2(); Vecteur evalyp(double x, Vecteur y); // integration par Runge-Kutta d'ordre 4 Vecteur rk4(double x0,Vecteur y0, double xf); void ecrirexyz(const char * nom_fichier); void ecrirexyzROOT(const char * nom_fichier); }; #endif
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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 #include "EDO2.h" #include <fstream> #include <iomanip> #include <iostream> using namespace std; #include "TGraph.h" #include "TFile.h" #include "TString.h" // On veux résoudre y"=-ay, on pose y1=y et y2=y' => y1'=y2 y2'=-ay1 EDO2::EDO2(int Ni, long double Ei) : N(Ni), E(Ei) { yf.redim(2); // On alloue les tableaux servant a stocker x, y, y', z, z' pour chaque point. xtab=new double [N+1]; // [0..N] => N+1 ytab=new Vecteur [N+1]; yptab=new Vecteur [N+1]; for (int i=0;i<N+1;i++) { ytab[i].redim(2); yptab[i].redim(2); } } EDO2::~EDO2() { delete[] xtab; delete[] ytab; delete[] yptab; } void EDO2::stocke(int i,double x, Vecteur y, Vecteur yp) { xtab[i]=x; ytab[i]=y; yptab[i]=yp; } void EDO2::ecrirexyz(const char * nom_fichier) { ofstream fichier( nom_fichier ); // On ecrit tout en 5 colonnes dans un fichier for (int i = 0 ; i < N+1 ; i++) { fichier << setw(20) << setprecision(12) << xtab[i]; for ( int d=1; d<=1; d++){ fichier << setw(20) << setprecision(12) << ytab[i](d)<< setw(20) << setprecision(12) << yptab[i](d); } fichier << endl;} fichier.close(); } void EDO2::ecrirexyzROOT(const char * nom_fichier) { TFile fichier( nom_fichier , "RECREATE" ); TGraph* gr=new TGraph [2]; TGraph* grp=new TGraph [2]; for (int d=0; d<1; d++) { TString s("Fonction "); TString name("F"); gr[d].Set(N+1); // mets le nombre de points dans les graphes s+=(d+1);name+=(d+1); gr[d].SetNameTitle(name,s); grp[d].Set(N+1); // mets le nombre de points dans les graphes des derivees s+=" derivee"; name+="p"; grp[d].SetNameTitle(name,s); for (int i = 0 ; i < N+1 ; i++) { gr[d].SetPoint(i,xtab[i],ytab[i](d+1)); grp[d].SetPoint(i,xtab[i],yptab[i](d+1)); } gr[d].Write(); grp[d].Write(); } delete [] gr; delete [] grp; fichier.Close(); } Vecteur EDO2::evalyp(double x, Vecteur y) { long double K, Ki=-((V-E)*2*m)/(hb*hb), Ke=-((-E)*2*m)/(hb*hb); //a=2m/h if (x>=0 && x<=L) K=Ki; else K=Ke; Vecteur v(2); v(1)=y(2); v(2)=-K*y(1); return v; } Vecteur EDO2::rk4(double x0, Vecteur y0,double xf) {double x; Vecteur y(2); Vecteur yp(2); Vecteur k1(2), k2(2), k3(2), k4(2); y=y0; double h = (xf - x0)/N; for (int i = 0 ; i<N ; i++) { x = x0 + i*h; yp= evalyp(x, y); stocke(i,x,y,yp); k1 = h * yp; k2 = h * evalyp(x+h/2, y+k1/2 ); k3 = h * evalyp(x+h/2, y+k2/2); k4 = h * evalyp(x+h, y+k3 ); y=y+(k1+2*k2+2*k3+k4)/6; } yf=y; stocke(N,xf,y,evalyp(xf,y)); return yf; }
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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 #include <iostream> #include <iomanip> #include <cmath> #include <vector> #include "EDO2.h" using namespace std; #include "TGraph.h" #include "TFile.h" #include <cmath> int main() { cout << endl << "EDO2" << endl; cout << endl << "N=1000" << endl; EDO2 T1(1000, 1.1e-19); Vecteur yres(2); double y00[2]={0,cos(T1.k1*100e-9)}; Vecteur y0(2,y00); yres=T1.rk4(100e-9,y0,-100e-9); T1.ecrirexyz("Diffusion.dat"); T1.ecrirexyzROOT("Diffusion.root"); //k1=k3=sqrt((2*m*E)/(h*h)) //k2=sqrt((2*m*(E-V))/(h*h)) // 1/T + R/T = phi0 // (1/T - R/T)k1 = phi0' double psi0, psi0p; psi0=T1.rk4(100e-9,y0,0)(1); psi0p=T1.rk4(100e-9,y0,0)(2); cout<<"psi0="<<psi0<<endl; cout<<"psi0p="<<psi0p<<endl; double A,B,R,T; A=0.5*(psi0+(psi0p/T1.k1)); B=psi0-A; cout<<"A="<<A<<endl; cout<<"B="<<B<<endl; T=1/A; R=B*T; cout<<"T="<<T<<endl; cout<<"R="<<R<<endl; cout<<"R+T="<<R+T<<endl; cout<<"k(A-B)="<<T1.k1*((1/T)-(R/T))<<endl; return 0;}
Partager