Table des matières PDFPython

Équation de Schrödinger à une dimension

1. Introduction

Ce document montre comment faire la résolution numérique de l'équation de Schrödinger dépendante du temps, en suivant une méthode similaire à celle utilisée pour l'équation de diffusion. On verra comment appliquer la méthode à un paquet d'onde diffusé par un potentiel.

Voir aussi la simulation Équation de Schrödinger à une dimension.

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 à la diffusion de particules dans un potentiel V(x). La résolution est limitée à l'intervalle fini [0,L]. Sur les bords de cet intervalle, la fonction d'onde est supposée nulle. La condition de normalisation s'écrit :

La fonction d'onde initiale est connue. Elle prendra la forme d'un paquet d'onde localisée à l'intérieur de l'intervalle [0,L].

On cherche à résoudre numériquement l'équation par la méthode des différences finies. La méthode utilisée est similaire à celle développée dans Équation de diffusion à une dimension.

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

et l'échelle de longueur :

En utilisant un temps et un x réduits par ces deux échelles, on obtient l'équation :

On utilise l'opérateur hamiltonien :

La solution de l'équation de Schrödinger peut s'écrire formellement à l'aide de l'exponentielle de l'opérateur hamiltonien ([1]) :

La conservation de la condition au cours du temps vient du caractère unitaire de l'opérateur .

3. Méthode numérique

3.a. Discrétisation

Soit N le nombre de points dans l'intervalle [0,L]. On définit le pas d'espace :

Soit Δt le pas de temps. On pose :

n est un entier positif ou nul repésentant le temps.

La dérivée temporelle est discrétisée de la manière suivante :

La dérivée seconde par rapport à x est discrétisée par :

Le schéma explicite qui en découle est :

Comme pour l'équation de diffusion, ce schéma explicite a une condition de stabilité très contraignante. On lui préfère donc un schéma implicite.

Pour alléger les notations, on introduit l'opérateur hamiltonien discret défini par :

3.b. Schéma implicite unitaire

En suivant la démarche adoptée pour l'équation de diffusion à une dimension, on est amené à considérer le schéma implicite suivant :

ou encore :

Ce schéma conduit à la résolution, pour chaque pas de temps, d'un système linéaire tridiagonal, pour lequel il existe une méthode directe efficace. Il a toutefois l'inconvénient de ne pas conserver la norme de la fonction d'onde.

Pour un intervalle de temps, on a :

Le schéma défini ci-dessus revient à utiliser l'approximation au premier ordre suivante :

Cette transformation n'est pas unitaire, c'est pourquoi la norme de la fonction d'onde n'est pas conservée. On la remplace par l'approximation de Cayley ([2]), qui est unitaire :

Le schéma numérique s'écrit alors :

On obtient ainsi le schéma implicite suivant :

avec :

On le met sous la forme générale suivante :

Pour chaque pas de temps, le système tridiagonal est résolu par l'algorithme d'élimination décrit dans Équation de diffusion à une dimension.

Les conditions limites sont de type Dirichlet, puisqu'elles consistent à annuler la fonction d'onde sur les bords de l'intervalle :

3.c. Paquet d'onde

Pour obtenir l'évolution d'un paquet d'onde, on considère comme condition initiale un paquet d'onde gaussien défini par :

Le paquet d'onde est initialement en x=x0 et se déplace dans le sens de l'axe x. La densité de propabilité (non normée) associée à ce paquet est :

La variance de cette densité de probabilité est σ0. Le spectre en nombre d'onde de ce paquet a une largeur donnée par :

Considérons l'évolution du paquet d'onde en l'absence de potentiel (particule libre). La largeur Δk reste constante au cours du temps. Notons σ la largeur du paquet à l'instant t. L'inégalité d'Heisenberg s'écrit :

Pour ce paquet gaussien, l'égalité est réalisée à l'instant t=0. On peut calculer précisément la largeur du paquet à l'instant t ([3]) :

Cette relation montre l'étalement du paquet d'onde lorsqu'on s'éloigne de l'origine du temps.

La vitesse du paquet d'onde (la vitesse de groupe) est :

Soit T la durée de la simulation. La distance parcourue par le paquet pendant cette durée doit être de l'ordre de la largeur L. Avec les unités utilisées (m=1/2 et ℏ=1), la condition s'écrit :

Si l'on souhaite que l'étalement du paquet d'onde soit négligeable pendant cette durée, il faut respecter la condition suivante :

qui s'écrit aussi :

Par exemple, pour L=1 et σ0=0.05, il faut que k0>200.

Voyons comment choisir l'échantillonnage spatial et temporel. Le pas d'espace doit vérifier :

La condition de Nyquist-Shannon impose une valeur maximale égale à π, mais il faut prévoir une marge à cause de la dispersion Δk. L'évolution dans le temps de la fonction d'onde se fait à la pulsation moyenne :

Avec les unités utilisées, le pas de temps doit donc vérifier :

Pour fixer les pas de temps et d'espace, on peut adopter la relation suivante :

Il suffit alors de vérifier la relation :

Par exemple, pour k0=200 et L=1, il faut N=2000 points. Le nombre de pas de temps nécessaire est :

Avec les valeurs précédentes P=5000.

On voit donc que la simulation d'un paquet d'onde présentant un faible étalement est très exigeante en terme de quantité de calculs.

4. Implémentation python

import numpy
import math
import cmath


class Schrodinger1D:
    def __init__(self,L,N):
        self.L = L
        self.N = N
        self.dx = L/(N-1)
        self.b = numpy.zeros(N,dtype=numpy.complex)
        self.e = numpy.zeros(N,dtype=numpy.complex)
        self.x = numpy.zeros(N,dtype=numpy.complex)
        self.beta = numpy.zeros(N,dtype=numpy.complex)
        self.gamma = numpy.zeros(N,dtype=numpy.complex)
        
    def config(self,dt,V):
        N = self.N
        self.dt = dt
        dx2 = self.dx*self.dx
        alpha = 2*dx2/dt
        for k in range(1,self.N-1):
            self.b[k] = 1j*alpha-dx2*V[k]-2.0
            self.e[k] = 1j*alpha+dx2*V[k]+2.0
        self.b[0] = 1.0
        self.b[N-1] = 1.0
        self.beta[0] = self.b[0]
        self.gamma[0] = 1.0/self.beta[0]
        for k in range(1,N):
            self.beta[k] = self.b[k]-self.gamma[k-1]
            if self.beta[k]==0:
                raise Exception("Impossible de resoudre le systeme ligne "+str(k))
            self.gamma[k] = 1.0/self.beta[k]
            
    def iterations(self,psi0,ti,tf):
        if psi0.size!=self.N:
            raise Exception("Taille de U incorrecte")
        psi = psi0.copy()
        t = ti
        while t<tf:
            self.x[0] = (self.e[0]*psi[0]-psi[1])/self.beta[0]
            for k in range(1,self.N-1):
                self.x[k] = (-psi[k-1]+self.e[k]*psi[k]-psi[k+1]-self.x[k-1])/self.beta[k]
            k = self.N-1
            self.x[k] = (-psi[k-1]+self.e[k]*psi[k]-self.x[k-1])/self.beta[k]
            psi[self.N-1] = self.x[self.N-1]
            for k in range(self.N-2,-1,-1):
                psi[k] = self.x[k]-self.gamma[k]*psi[k+1]
            t += self.dt
        return [psi,t]
            
    def V_marche(self,V0):
        V = numpy.zeros(self.N)
        P = int(self.N/2)
        V[P:self.N-1] = V0
        return V
        
    def V_barriere(self,V0,largeur):
        q = int(largeur/self.dx/2)
        V = numpy.zeros(self.N)
        P = int(self.N/2)
        V[P-q:P+q] = V0
        return V
     
    def paquet(self,x0,k0,sigma0):
        psi = numpy.zeros(self.N,dtype=numpy.complex)
        for k in range(1,self.N-1):
            x = self.dx*k
            psi[k] = cmath.exp(1j*k0*x)*math.exp(-(x-x0)*(x-x0)/(4*sigma0**2))
        return psi
            

5. Exemples

5.a. Diffusion de particules sur une marche de potentiel

import numpy
from Schrodinger1D import Schrodinger1D
from matplotlib.pyplot import *

N=2000
L=1.0
dx=L/(N-1)
x = numpy.arange(N)*dx
schrodinger = Schrodinger1D(L,N)
V = schrodinger.V_marche(1.0)
l0 = 2e-2
k0 = numpy.pi*2/l0
E=k0**2
x0 = 0.25
sigma0 = 0.04
psi = schrodinger.paquet(x0,k0,sigma0)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()


                
figAfigA.pdf

On trace le module de la fonction d'onde (densité de probabilité) pour différents instants.

t=0.0
dt=2*dx*dx
v=2*k0
T = 1.0/v
V0 = E
schrodinger.config(dt,V*V0)
[psi,t] = schrodinger.iterations(psi,t,T*0.2)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figBfigB.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.3)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figCfigC.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.5)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figDfigD.pdf

5.b. Diffusion de particules sur une barrière de potentiel

N=2000
L=1.0
dx=L/(N-1)
x = numpy.arange(N)*dx
schrodinger = Schrodinger1D(L,N)
V = schrodinger.V_barriere(1.0,0.05)
l0 = 2e-2
k0 = numpy.pi*2/l0
E=k0*k0
x0 = 0.25
sigma0 = 0.04
psi = schrodinger.paquet(x0,k0,sigma0)
                

On trace le module de la fonction d'onde (densité de probabilité) pour différents instants.

t=0.0
dt=2*dx*dx
v=2*k0
T = 1.0/v
V0=E
schrodinger.config(dt,V*V0)
[psi,t] = schrodinger.iterations(psi,t,T*0.2)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figEfigE.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.3)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figFfigF.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.5)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figGfigG.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.55)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figHfigH.pdf
Références
[1]  J.J. Sakurai,  Modern quantum mechanics,  (Addison-Wesley, 1985)
[2]  A. Goldberg, H.M. Schey, J.L. Schwartz,  Computer-generated motion pictures of one-dimensional quantum-mechanical transmission and reflection phenomena,  (Am. J. Phys, 35(3), 1967)
[3]  C. Cohen-Tannoudji, B. Diu, F. Laloë,  Mécanique quantique,  (Hermann, 1980)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.