Table des matières PDF

Équation de Schrödinger à deux dimensions

1. Introduction

2. Équation de Schrödinger

On considère l'équation de Schrödinger dépendante du temps pour une fonction d'onde :

On s'intéresse au mouvement d'une particule dans un potentiel, limité par le domaine rectangulaire [0,Lx]x[0,Ly]. Sur les bords de ce domaine, la fonction d'onde est nulle. La condition de normalisation s'écrit :

La fonction d'onde initiale est connue, et prendra la forme d'un paquet d'onde localisé dans le domaine.

L'objectif est de résoudre numériquement l'équation par la méthode des différences finies, en suivant une méthode similaire à celle exposée dans Équation de diffusion à deux dimensions. L'équation de Schrödinger est en effet analogue à l'équation de diffusion en raison de la dérivée première par rapport au temps. Elle a des solutions sous forme d'onde progressives en raison du coefficient imaginaire .

Soit V0 l'échelle du potentiel. La fonction V(x,y) utilisée sera sans dimensions. On introduit l'échelle de temps :

et l'échelle de longueur :

En utilisant un temps et des coordonnées (x,y) réduits par ces deux échelles, on obtient l'équation :

On utilise l'opérateur hamiltonien :

3. Méthode numérique

3.a. Discrétisation

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

L'indice j varie de 0 à Nx-1, l'indice m de 0 à Ny-1.

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

n est un entier positif ou nul. On pose :

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 :

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

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 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 :

On a bien sûr :

Le premier demi-pas s'écrit :

Le second demi-pas s'écrit :

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

Le premier demi-pas s'écrit :

Le second demi-pas s'écrit :

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 :

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

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

et tous les autres coefficients dont nuls.

3.c. Algorithme

La maillage est constitué de matrices à Nx colonnes et Ny 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 à Nx :

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

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 à Nx:

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 Nx à 1 :

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 à Ny :

puis pour j variant de Ny à 1 :

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 Nx=23 et Ny=23 (les points sont les nœuds de la grille).

matrice.svgFigure pleine page

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

Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.