Équation de Schrödinger à deux dimensions

3.a. Discrétisation

Les nœuds du maillage sont définis par :

$$\begin{array}{ccc}\hfill & x_j=j\Delta x\hfill & \hfill \text{(7)}\\ \hfill & y_m=m\Delta y\hfill & \hfill \text{(8)}\end{array}$$

L'indice j varie de 0 à N_x-1, l'indice m de 0 à N_y-1.

Les valeurs sont recherchées aux instants définis par :

$$t_n=n\Delta t\text{(9)}$$

où n est un entier positif ou nul. On pose :

$$\begin{array}{ccc}\hfill & {\psi }_{j,m}^n=\psi (j\Delta x,m\Delta y,n\Delta t)\hfill & \hfill \text{(10)}\\ \hfill & V_{j,m}=V(j\Delta x,m\Delta y)\hfill & \hfill \text{(11)}\end{array}$$

Dans cette notation, l'exposant sert donc à représenter le temps.

La discrétisation des dérivées secondes par la méthode des différences finies nous amène à définir un opérateur hamiltonien :

$$H{\psi }_{j,m}^n=-\frac{{\psi }_{j+1,m}^n-2{\psi }_{j,m}^n+{\psi }_{j-1,m}^n}{{\left(\Delta x\right)}^2}-\frac{{\psi }_{j,m+1}^n-2{\psi }_{j,m}^n+{\psi }_{j,m-1}^n}{{\left(\Delta y\right)}^2}+V_{j,m}{\psi }_{j,m}^n\text{(12)}$$

Comme il est expliqué dans Équation de Schrödinger à une dimension, un schéma numérique qui conserve la condition de normalisation s'écrit :

$$(1+i\frac{1}{2}\Delta tH){\psi }_{j,m}^{n+1}=(1-i\frac{1}{2}\Delta tH){\psi }_{j,m}^n\text{(13)}$$

Le schéma numérique obtenu est implicite.

3.b. Schéma numérique à directions alternées

Il s'agit de modifier le schéma implicite $\mathrm{(13)}$ pour obtenir un schéma implicite dans la direction x, alterné avec un schéma implicite dans la direction y, comme exposé dans Équation de diffusion à deux dimensions.

Le schéma implicite à directions alternées consiste à décomposer le pas de temps (n,n+1) en deux sous-pas. Pour le premier, on écrit un schéma implicite pour la dérivée par rapport à x mais explicite pour la dérivée par rapport à y. Pour le second, on écrit un schéma implicite pour la dérivée par rapport à y mais explicite pour la dérivée par rapport à x. On choisit de prendre en compte le potentiel seulement dans le second demi-pas.

Pour alléger les notations, on définit les opérateurs de dérivation suivants :

$$\begin{array}{ccc}\hfill & D_x{\psi }_{j,m}^n=\frac{{\psi }_{j+1,m}^n-2{\psi }_{j,m}^n+{\psi }_{j-1,m}^n}{{\left(\Delta x\right)}^2}\hfill & \hfill \text{(14)}\\ \hfill & D_y{\psi }_{j,m}^n=\frac{{\psi }_{j,m+1}^n-2{\psi }_{j,m}^n+{\psi }_{j,m-1}^n}{{\left(\Delta y\right)}^2}\hfill & \hfill \text{(15)}\end{array}$$

On a bien sûr :

$$H=-D_x-D_y+V_{j,m}\text{(16)}$$

Le premier demi-pas s'écrit :

$$(1-i\frac{\Delta t}{2}{D}_x+i\frac{\Delta t}{2}\frac{1}{2}{V}_{j,m}){\psi }_{j,m}^{n+\frac{1}{2}}=(1+i\frac{\Delta t}{2}{D}_y-i\frac{\Delta t}{2}\frac{1}{2}{V}_{j,m}){\psi }_{j,m}^n\text{(17)}$$

Le second demi-pas s'écrit :

$$(1-i\frac{\Delta t}{2}{D}_y+i\frac{\Delta t}{2}\frac{1}{2}{V}_{j,m}){\psi }_{j,m}^{n+1}=(1+i\frac{\Delta t}{2}{D}_x-i\frac{\Delta t}{2}\frac{1}{2}{V}_{j,m}){\psi }_{j,m}^{n+\frac{1}{2}}\text{(18)}$$

Pour simplifier, on pose Δx=Δy et on définit :

$$\alpha =\frac{2{\left(\Delta x\right)}^2}{\Delta t}\text{(19)}$$

Le premier demi-pas s'écrit :

$${\psi }_{j,m}^{n+\frac{1}{2}}-\frac{i}{\alpha }({\psi }_{j+1,m}^{n+\frac{1}{2}}-2{\psi }_{j,m}^{n+\frac{1}{2}}+{\psi }_{j-1,m}^{n+\frac{1}{2}})+\frac{i}{\alpha }{\left(\Delta x\right)}^2\frac{1}{2}{V}_{j,m}{\psi }_{j,m}^{n+\frac{1}{2}}={\psi }_{j,m}^n+\frac{i}{\alpha }({\psi }_{j,m+1}^n-2{\psi }_{j,m}^n+{\psi }_{j,m-1}^n)-\frac{i}{\alpha }{\left(\Delta x\right)}^2\frac{1}{2}{V}_{j,m}{\psi }_{j,m}^n\text{(20)}$$

Le second demi-pas s'écrit :

$${\psi }_{j,m}^{n+1}-\frac{i}{\alpha }({\psi }_{j,m+1}^{n+1}-2{\psi }_{j,m}^{n+1}+{\psi }_{j,m-1}^{n+1})+\frac{i}{\alpha }{\left(\Delta x\right)}^2\frac{1}{2}{V}_{j,m}{\psi }_{j,m}^{n+1}={\psi }_{j,m}^{n+\frac{1}{2}}+\frac{i}{\alpha }({\psi }_{j+1,m}^{n+\frac{1}{2}}-2{\psi }_{j,m}^{n+\frac{1}{2}}+{\psi }_{j-1,m}^{n+\frac{1}{2}})-\frac{i}{\alpha }{\left(\Delta x\right)}^2\frac{1}{2}{V}_{j,m}{\psi }_{j,m}^{n+\frac{1}{2}}\text{(21)}$$

Chacune de ces deux étapes consiste à résoudre un système tridiagonal, ce qui se fait par élimination directe en utilisant l'algorithme décrit dans équation de diffusion à une dimension.

On écrira les deux systèmes à résoudre sous la forme générale suivante :

$$\begin{array}{ccc}\hfill & A_{j,m}{\psi }_{j-1,m}^{n+1/2}+B_{j,m}{\psi }_{j,m}^{n+1/2}+C_{j,m}{\psi }_{j+1,m}^{n+1/2}=D_{j,m}{\psi }_{j,m-1}^n+E_{j,m}{\psi }_{j,m}^n+F_{j,m}{\psi }_{j,m+1}^n+G_{j,m}\hfill & \hfill \text{(22)}\\ \hfill & A{\text{'}}_{j,m}{\psi }_{j,m-1}^{n+1}+B{\text{'}}_{j,m}{\psi }_{j,m}^{n+1}+C{\text{'}}_{j,m}{\psi }_{j,m+1}^{n+1}=D{\text{'}}_{j,m}{\psi }_{j-1,m}^{n+1/2}+E{\text{'}}_{j,m}{\psi }_{j,m}^{n+1/2}+F{\text{'}}_{j,m}{\psi }_{j+1,m}^{n+1/2}+G{\text{'}}_{j,m}\hfill & \hfill \text{(23)}\end{array}$$

Pour un point non situé sur un bord, les coefficients sont :

$$\begin{array}{ccc}\hfill & A_{j,m}=A{\text{'}}_{j,m}=1\hfill & \hfill \text{(24)}\\ \hfill & B_{j,m}=B{\text{'}}_{j,m}=i\alpha -2-{\left(\Delta x\right)}^2\frac{1}{2}{V}_{j,m}\hfill & \hfill \text{(25)}\\ \hfill & C_{j,m}=C{\text{'}}_{j,m}=1\hfill & \hfill \text{(26)}\\ \hfill & D_{j,m}=D{\text{'}}_{j,m}=-1\hfill & \hfill \text{(27)}\\ \hfill & E_{j,m}=E{\text{'}}_{j,m}=i\alpha +2+{\left(\Delta x\right)}^2\frac{1}{2}{V}_{j,m}\hfill & \hfill \text{(28)}\\ \hfill & F_{j,m}=F{\text{'}}_{j,m}=-1\hfill & \hfill \text{(29)}\\ \hfill & G_{j,m}=G{\text{'}}_{j,m}=0\hfill & \hfill \text{(30)}\end{array}$$

Sur les bords du domaine, ou sur tout point où une valeur nulle est attribuée à la fonction d'onde, on a :

$$B_{j,m}=B{\text{'}}_{j,m}=1\text{(31)}$$

et tous les autres coefficients dont nuls.

3.c. Algorithme

La maillage est constitué de matrices à N_x colonnes et N_y lignes. La définition des pas de temps et d'espace, des conditions limites et des sources permet de construire les matrices A,B,C,D,E,F,G,A',B',C',D',E',F',G'.

Il faut aussi construire les matrices β,γ pour la résolution des systèmes sur les lignes, et les matrices β',γ' pour la résolution des systèmes sur les colonnes. Pour chaque ligne j, la récurrence à appliquer pour i variant de 1 à N_x :

$$\begin{array}{ccc}\hfill & {\beta }_{1,m}=B_{1,m}\hfill & \hfill \text{(32)}\\ \hfill & {\gamma }_{1,m}=\frac{C_{1,m}}{{\beta }_{1,m}}\hfill & \hfill \text{(33)}\\ \hfill & {\beta }_{j,m}=B_{j,m}-A_{j,m}{\gamma }_{i-1,m}\hfill & \hfill \text{(34)}\\ \hfill & {\gamma }_{j,m}=\frac{C_{j,m}}{{\beta }_{j,m}}\hfill & \hfill \text{(35)}\end{array}$$

Pour chaque colonne i, la récurrence à appliquer pour j variant de 1 à N_y :

$$\begin{array}{ccc}\hfill & \beta {\text{'}}_{j,1}=B{\text{'}}_{j,1}\hfill & \hfill \text{(36)}\\ \hfill & \gamma {\text{'}}_{j,1}=\frac{C{\text{'}}_{j,1}}{\beta {\text{'}}_{j,1}}\hfill & \hfill \text{(37)}\\ \hfill & \beta {\text{'}}_{j,m}=B{\text{'}}_{j,m}-A{\text{'}}_{j,m}\gamma {\text{'}}_{j,m-1}\hfill & \hfill \text{(38)}\\ \hfill & \gamma {\text{'}}_{j,m}=\frac{C{\text{'}}_{j,m}}{\beta {\text{'}}_{j,m}}\hfill & \hfill \text{(39)}\end{array}$$

Pour chaque pas de temps n, on doit effectuer une résolution tridiagonale pour chaque ligne, puis pour chaque colonne.

Pour chaque ligne j on commence par calculer la suite suivante pour i variant de 1 à N_x:

$$\begin{array}{ccc}\hfill & x_1=\frac{1}{{\beta }_{1,m}}(D_{1,m}{\psi }_{1,m-1}^n+E_{1,m}{\psi }_{1,m}^n+F_{1,m}{\psi }_{1,m+1}^n+G_{1,m})\hfill & \hfill \text{(40)}\\ \hfill & x_i=\frac{1}{{\beta }_{j,m}}(D_{j,m}{\psi }_{j,m-1}^n+E_{j,m}{\psi }_{j,m}^n+F_{j,m}{\psi }_{j,m+1}^n+G_{j,m}-A_{j,m}{x}_{i-1})\hfill & \hfill \text{(41)}\end{array}$$

Pour la même ligne j, on calcule les nouvelles valeurs à l'instant n+1/2 par récurrence, en faisant varier de i de N_x à 1 :

$$\begin{array}{ccc}\hfill & {\psi }_{N_x,m}^{n+\frac{1}{2}}=x_{N_x}\hfill & \hfill \text{(42)}\\ \hfill & {\psi }_{j,m}^{n+\frac{1}{2}}=x_i-{\gamma }_{j,m}{\psi }_{i+1,m}^{n+\frac{1}{2}}\hfill & \hfill \text{(43)}\end{array}$$

Les valeurs à l'instant n+1/2 sont stockées dans la matrice ψ. On procède alors au parcours des colonnes. Pour chaque colonne i, on calcule la suite suivante pour j variant de 1 à N_y :

$$\begin{array}{ccc}\hfill & x_1=\frac{1}{\beta {\text{'}}_{j,1}}(D{\text{'}}_{j,1}{\psi }_{i-1,1}^{n+1/2}+E{\text{'}}_{j,1}{\psi }_{j,1}^{n+1/2}+F{\text{'}}_{j,1}{\psi }_{i+1,1}^{n+1/2}+G{\text{'}}_{j,1})\hfill & \hfill \text{(44)}\\ \hfill & x_i=\frac{1}{\beta {\text{'}}_{j,m}}(D{\text{'}}_{j,m}{\psi }_{i-1,m}^{n+1/2}+E{\text{'}}_{j,m}{\psi }_{j,m}^{n+1/2}+F{\text{'}}_{j,m}{\psi }_{i+1,m}^{n+1/2}+G{\text{'}}_{j,m}-A{\text{'}}_{j,m}{x}_{i-1})\hfill & \hfill \text{(45)}\end{array}$$

puis pour j variant de N_y à 1 :

$$\begin{array}{ccc}\hfill & {\psi }_{j,N_y}^n=x_{N_y}\hfill & \hfill \text{(46)}\\ \hfill & {\psi }_{j,m}^n=x_i-\gamma {\text{'}}_{j,m}{\psi }_{j,m+1}^n\hfill & \hfill \text{(47)}\end{array}$$

Lorsqu'on traite les lignes, il faut écrire dans une seconde matrice ψ, pour ne pas écraser les valeurs de l'instant n. Il en est de même lorsqu'on traite les colonnes.

Le nombre de matrices utilisées est 20. Pour un maillage 1024x1024 en flottants 32 bits, la mémoire utilisée est de 80 méga-octets.

Les lignes peuvent être parcourues en parallèle (idem pour les colonnes). On peut donc avec profit utiliser un processeur graphique pour faire ces calculs (via l'API CUDA ou OpenCL). Pour le traitement des lignes, chaque unité de calcul traite une ligne. Les opérations effectuées sur les différentes lignes sont rigoureusement identiques, ce qui est le cas idéal pour une implémentation sur GPU. Bien entendu, toutes les matrices doivent être transférées sur la mémoire de la carte graphique avant le calcul.

On remarque que le traitement des bords (qui doit être identique au traitement des autres points) nécessite un accès à des éléments de ψ hors domaine. Pour résoudre ce problème, on utilise des matrices avec des rangées fictives sur les bords. Le nombre d'éléments dans chaque direction est une puissance de 2 (plus efficace pour l'implémentation GPU). La figure suivante montre un exemple pour N_x=2³ et N_y=2³ (les points sont les nœuds de la grille).

Figure pleine page

Sur cet exemple, les matrices sont 10x10 mais le maillage effectif est 8x8.