Chapitre3 Équations aux dérivées partielles et

MASTER PHYSIQUE - Toutes options confondues

fonctions de Green
Samir Kenouche - Département des Sciences de la Matière - UMKB
Méthodes Mathématiques et Algorithmes pour la Physique
Résumé
Ce chapitre débute par un rappel succinct de quelques notions élémentaires sur les équations différentielles.
Les méthodes numériques d’Euler, de Heun, de Crank-Nicolson et de Runge-Kutta ont été étudiées en détail
en deuxième années L2, module : Méthodes numériques et programmation. A cet effet, je renvoie les lecteurs
intéressés à mon site web personnel. Ainsi, il a été jugé inutile de les rappeler ici, en revanche, l’accent a été
mis particulièrement sur les résolutions analytique et numérique des équations aux dérivées partielles. A cet
égard, de nombreux problèmes issus de la physique ont été résolus, à l’instar de : l’équation de Schrödinger, la
corde élastique, diffusion de la chaleur, convection-diffusion de la chaleur et l’équation de la flexion. L’équation
de Schrödinger est l’équation maitresse de la mécanique quantique, décrivant les propriétés de la matière à
des échelles atomiques. Cette équation a été résolue en détail sans occultation des interprétations physiques
des différentes étapes. Les autres problèmes (linéaire et non-linéaire) ont été résolus numériquement, au
moyen de scripts Matlabr , en considérant la méthode des différences finies. Finalement, soulignons que les
fondements mathématiques des méthodes numériques sont délibérément occultés afin de privilégier d’avantage
les aspects opérationnels, qui sont les plus importants pour un non-mathématicien. L’objectif est de conférer
aux étudiants (es) issus des sciences expérimentales, un certain nombre d’outils mathématiques, afin qu’ils
puissent savoir se servir d’un ordinateur pour résoudre différents problèmes de la physique.

Table des matières

Cours complet est disponible sur mon site web : [Link]

I Introduction 58
I-A Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

II Équation aux dérivées partielles 60

II-A Résolution analytique de quelques EDPs . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
II-A1 Établissement de l’équation de Schrödinger . . . . . . . . . . . . . . . . . . . . . 61
II-A2 Équation dépendante du temps . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
II-B Solution exacte de l’équation de Schrödinger . . . . . . . . . . . . . . . . . . . . . . . . 66
II-B1 Résolution de la partie angulaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
II-B2 Résolution de la partie radiale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
II-C Résolution numérique de quelques EDPs . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
II-C1 En dimension 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
II-C2 Problème non-linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
II-C3 En dimension 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

III Fonctions de Green 93

III-A Fonction de Green de l’équation de Schrödinger . . . . . . . . . . . . . . . . . . . . . . 99

IV Annexe : Compléments sur les orbitales des atomes réels 106

IV-A Orbitales de type Slater . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
IV-B Orbitales Gaussienne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
S. Kenouche est docteur en Physique de l’Université des Sciences et Techniques de Montpellier (France) et docteur en Chimie
de l’Université de Béjaia (Algérie).
Site web : voir [Link]
Version corrigée, améliorée et actualisée le 10.10.2020.
 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

I. Introduction

U Ne équation différentielle est une relation fonctionnelle dont l’inconnue est une fonction y(t), avec t ∈
[a, b]. La forme générale d’une telle équation s’écrit :

f (y (n) , y (n−1) , y (n−2) , ..., y (1) , y, t) = φ(t) (1)

Avec, y (n) est la nième dérivée de la fonction y et ϕ(t) désigne le second membre de l’équation différentielle.
Dans le cas où φ(t) = 0, on dira que l’équation différentielle est homogène. L’existence d’une solution unique
de l’équation différentielle est tributaire de l’imposition de certaines conditions initiales sur y(t) et ses dérivées.
Dans l’équation ci-dessus, les conditions initiales sont les valeurs de y(a), y (1) (a), y (2) (a), ..., y (n) (a). Cependant,
il faut noter que très souvent la solution analytique n’existe pas, et on doit par conséquent approcher la solution
exacte y(t) par des méthodes numériques en calculant des approximations successives à chaque instant.

A. Problème de Cauchy
Le problème de Cauchy consiste à trouver une fonction continûment dérivable y : t ∈ R+ −→ y(t) ∈ R
vérifiant :
 0
 y (t) = f (t, y(t)),
 t>0
P: (2)
 y(0) = y

0

La première équation est une équation différentielle et la deuxième relation exprime une condition de Cauchy
ou condition initiale. Une solution y(t) au problème (2) est appelée intégrale de l’équation différentielle. En
effet, ce problème est équivalent à l’équation intégrale :
Z t
y(t) = y0 + f (u, y(u)) du (3)
t0

a) Définition: soit f : I × R 7−→ R une fonction donnée, s’il existe une constante L > 0 telle que :
Cours complet est disponible sur mon site web : [Link]

|f (t, u) − f (t, v)| ≤ L |u − v| (4)

∀ u, v ∈ R et ∀ t ∈ I alors f est dite lipschitzienne de rapport L sur I × R ou simplement L-lipschitzienne.

b) Théorème: si f est continue sur I × R et L-lipschitzienne par rapport à sa deuxième variable y(t)
alors le problème de Cauchy admet une solution unique sur I, ∀ u(0) ∈ R.

c) Démonstration: pour démontrer ce théorème nous considérons l’application ϕ : (R → Rn ) 7−→

(R → Rn ) qui est définie par :
Z t
ϕ(y)(t) = y0 + f (u, y(u)) du (5)
t0

De sorte que,
Z t0
ϕ(y)(t0 ) = y0 + f (u, y(u)) du = y0 (6)
t0
| {z }
=0

Ainsi ϕ(y)(t0 ) satisfait toujours la condition initiale. En outre, si elle admet un point fixe, ie :
Z t
ϕ(y) = y ⇒ y(t) = y0 + f (u, y(u)) du (7)
t0

dy(t) d
Z t 
0
⇒ =0+ f (u, y(u)) du ⇒ y (t) = f (t, y(t)) (8)
dt dt t0

Chapitre3 Équations aux dérivées partielles et fonctions de Green 58

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Donc un point fixe de l’application ϕ sera forcément une solution de l’équation différentielle avec la même
condition initiale. La démonstration du théorème revient désormais à prouver que l’application ϕ est contrac-
tante. Cela signifie qu’elle va contracter l’espace des fonctions de sorte à rapprocher toute fonction d’une
fonction qui est solution du problème. Ci-dessous la démonstration :

Z t Z t
(1) (1)
k ϕ (y1 )(t) − ϕ (y2 )(t) k = k y0 + f (u, y1 (u)) du − y0 − f (u, y2 (u)) du k inég. triang.
t0 t0
Z t
≤ k f (u, y1 (u)) − f (u, y2 (u)) k du avec f est L-lipschitzienne
t0
Z t
≤L k y1 (u) − y2 (u) k du
t0 | {z }
≤ ky1 −y2 k∞

≤ L (t − t0 ) k y1 − y2 k∞
En itérant une deuxième fois :

Z t Z t
(2) (2)
k ϕ (y1 )(t) − ϕ (y2 )(t) k = k y0 + f (u, ϕ(y1 )(u)) du − y0 − f (u, ϕ(y2 )(u)) du k
t0 t0
Z t
≤ k f (u, ϕ(y1 )(u)) − f (u, ϕ(y2 )(u)) k du
t0
Z t
≤ k ϕ(y1 )(u) − ϕ(y2 )(u) k du
t0
Z t
≤ L2 k y1 − y2 k∞ du
t0
(t − t0 )2
≤ L2k y1 − y 2 k ∞
2
Cours complet est disponible sur mon site web : [Link]

Par récurrence à la nième itération, nous obtenons :

(t − t0 )n
k ϕ(n) (y1 )(t) − ϕ(n) (y2 )(t) k ≤ Ln k y1 − y2 k ∞ (9)
n
(t − t0 )n
Le terme Ln tend vers zéro quand n tend vers l’infinie ⇒ l’application ϕ est contractante. Ce
n
théorème garantit l’unicité des solutions des équations différentielles pour une condition initiale donnée. Au-
trement dit, à deux conditions initiales différentes correspondent deux solutions différentes. Ceci est d’une
importance majeure pour pouvoir prédire l’état d’un système à un instant ultérieur. Comme par exemple la
prédiction de la trajectoire d’une particule à partir de l’instant courant.

Nous rappelons que pour des fonctions continues f : D ⊂ R 7−→ R ∀ p ≥ 1, les normes p sont définies par :
Z 1/p
p
k f kp = |f (t)| dt (10)
t∈D

En l’occurrence la norme Euclidienne p = 2 est définie par :

Z 1/2
2
k f k2 = |f (t)| dt (11)
t∈D

La norme infinie s’écrit :

k f k∞ = sup |f (t)| = lim |f (t)|p (12)

t∈D p−→+∞

Chapitre3 Équations aux dérivées partielles et fonctions de Green 59

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

II. Équation aux dérivées partielles
Une équation aux dérivées partielles (EDP en abrégé) exprime une relation fonctionnelle entre les variables
indépendantes X = (x, x, z, · · · ), la fonction inconnue u(x, y, z, · · · ) ∈ Ω ⊂ R et ses dérivées partielles. La
formulation générale d’une telle relation s’écrit suivant :
n
X ∂k
ak (x, y, z, · · · ) · Dk u(x, y, z, · · · ) = φ(x, y, z, · · · ), Dk ≡ (13)
k=0
∂X k
Ou de façon équivalente,
h i
F Dn u(x, y, z, · · · ), Dn−1 u(x, y, z, · · · ), · · · , D1 u(x, y, z, · · · ), u(x, y, z, · · · ), x, y, z, · · · = φ(x, y, z, · · · ) (14)
Où F exprime une relation fonctionnelle entre les variables indépendantes X = (x, x, z, · · · ), la fonction
inconnue u(x, y, z, · · · ) et ses dérivées partielles. On appelle ordre d’une EDP l’ordre de la plus grande dérivée
présente dans l’équation. On appelle problème aux limites, une EDP munie de conditions aux limites 1 sur
la totalité de la frontière du domaine sur lequel elle est posée. En se limitant aux EDPs de l’ordre deux
(k = 2), l’EDP exprimée par (13) est dite linéaire car les coefficients ak (x, y) 6= f (u) et également la source
φ(x, y) 6= f (u). Cette linéarité peut s’exprimer sous la forme : si u1 (x, y) et u2 (x, y) sont solutions de (13)
alors :

F [a1 u1 (x, y) + a2 u2 (x, y)] = a1 F [u1 (x, y)] + a2 F [u2 (x, y)] (15)
De la même façon, ∀λ ∈ R,

F [λ u(x, y)] = λ F [u(x, y)] (16)

Une EDP est dite semi-linéaire si elle est écrite sous la forme :
X
ak (x, y) · D(2) u(x, y) = φ(D(1) u(x, y), x, y) (17)
k≤2
Cours complet est disponible sur mon site web : [Link]

Une EDP est dite quasi-linéaire si elle est écrite sous la forme :
X
ak (D(1) u(x, y), u(x, y), x, y) · D(2) u(x, y) = φ(D(1) u(x, y), x, y) (18)
k≤2

Elle est linéaire seulement par rapport aux dérivées partielles de u(x, y) d’ordre deux. Une EDP linéaire
d’ordre deux est explicitée sous la forme :
∂ 2 u(x, y) ∂ 2 u(x, y) ∂ 2 u(x, y)
A(x, y) + B(x, y) + C(x, y) + (19)
∂x2 ∂x y ∂y 2
∂u(x, y) ∂u(x, y)
D(x, y) + E(x, y) + F (x, y)u(x, y) = φ(x, y)
∂x ∂y
Les coefficients peuvent êtres égaux à des constantes et A, B, C sont obligatoirement différents de zéro,
sinon on retombe sur une EDP d’ordre un. Notons que cette EDP est inhomogène car le second membre est
non nul, elle est dite homogène dans le cas contraire. Les EDPs linéaires d’ordre deux sont classées selon la
valeur du discriminant ∆, soit :


 ∆ = 0 ⇒ EDP parabolique




∆ = B(x, y)2 − 4 A(x, y) C(x, y) : ∆ > 0 ⇒ EDP hyperbolique (20)






∆ < 0 ⇒ EDP elliptique

1. Le lecteur est invité à distinguer les conditions dites de Dirichlet (ocndition imposée sur la valeur de la solution à la
frontière), de Neuman (condition imposée sur la valeur de la dérivée de la solution ) et de Dirichle-Neuman (condition mixte).
La précision globale du schéma dépend de la précision sur la discrétisation des conditions aux limites.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 60

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Il est possible de démontrer que cette classification est invariante par changement de bases dans le plan.
A titre informatif, les équations de Poisson et de Laplace sont elliptiques, celle de la chaleur est parabolique,
tandis que les équations des cordes vibrantes sont hyperboliques.

A. Résolution analytique de quelques EDPs

Dans ce qui suit, nous nous proposons de résoudre analytiquement l’équation des ondes ainsi que l’équation
maitresse de la théorie quantique, à savoir l’équation de Schrödinger. L’équation des ondes est résolue en
utilisant la méthode de séparation de variables (ou méthode de Fourier) et l’équation de Schrödinger est
résolue en se servant des séries entières.

1) Établissement de l’équation de Schrödinger: cette équation joue un rôle fondamental en mécanique

quantique car sa résolution permet la description des propriétés de la matière à des échelles atomiques. L’atome
d’hydrogène est l’un des rares systèmes réalistes (communément appelé système à deux corps) qui peuvent
réellement être résolus de manière analytique. La géométrie sphérique de ce système suggère l’utilisation
des coordonnées sphériques avec le noyau à l’origine. Erwin Schrödinger et Werner Heisenberg ont ébauché
séparément l’équation régissant la description et l’évolution des systèmes quantiques. Schrödinger a opté pour
un formalisme mathématique utilisant les équations aux dérivées partielles, alors que Heisenberg a utilisé un
formalisme matriciel. Bien que les deux approches se sont révélées mathématiquement équivalentes. La plupart
des ouvrages débutent par l’équation de Schrödinger, qui semble avoir une meilleure interprétation physique par
le biais de l’équation des ondes classique. En effet, l’équation de Schrödinger peut être vue comme une forme
de l’équation des ondes appliquée aux ondes de matière. L’équation des ondes classique unidimensionnelle est
donnée par :

∂ 2 u(x, t) 1 ∂ 2 u(x, t)
= (21)
∂ x2 v 2 ∂ t2
Avec les conditions aux limites :
Cours complet est disponible sur mon site web : [Link]

u(0, t) = 0 et u(l, t) = 0 (22)

Ces conditions stipulent que l’amplitude de vibration est nulle aux extrémités x = 0 et x = l. En séparant
les variables spatiale et temporelle :

u(x, t) = ψ(x)f (t)

Nous obtenons,

1 d2 ψ(x) 1 d2 f (t)
= (23)
ψ(x) dx2 v 2 f (t) dt2
Les deux termes sont égaux si et seulement s’ils valent la même constante, notée k. Cette dernière est
arbitraire dans le sens où son signe est inconnu. Elle peut être négative, positive ou même nulle :

1 d2 ψ(x)

=k

00
 
 ψ(x) − k ψ(x) = 0

 ψ(x) dx2

 
que nous écrirons sous la forme (24)


 1 d2 f (t)  f (t)00 − k v 2 f (t) = 0


 =k
v 2 f (t) dt2


En utilisant la méthode de séparation des variables, nous sommes passés d’une équation aux dérivées
partielles à deux équations différentielles ordinaires de second ordre sans second membre. La constante k

Chapitre3 Équations aux dérivées partielles et fonctions de Green 61

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

est appelée constante de séparation. Dans ce qui suit, nous résolvons l’équation différentielle dépendante de la
variable spatiale :
00
ψ(x) − k ψ(x) = 0 (25)
Afin d’atteindre cet objectif, nous envisageons trois cas de figure pour la constante de séparation k.

o Cas où k = 0

00 0
k = 0 ⇒ ψ(x) = 0 ⇒ ψ(x) = a1 ⇒ ψ(x) = a1 x + b1 (26)
Avec a1 et b1 sont des constantes d’intégration pouvant êtres déterminées tenant compte des conditions aux
limites (22) :
 
 u(0, t) = ψ(0) f (t) = 0
  ψ(0) = 0

f (t) 6= 0 ⇒ (27)
 u(l, t) = ψ(l) f (t) = 0
  ψ(l) = 0


A partir de (26) on obtient :

ψ(x = 0) = b1 = 0 et ψ(x = l) = a1 l + 0 = 0 ⇒ a1 = b1 = 0 ⇒ ψ(x) = 0

Autrement dit,

ψ(x) = 0 ⇒ u(x, t) = 0, ∀x ∈ [0 l]
Cette solution est mathématiquement juste mais physiquement inacceptable dans le sens où elle n’apporte
aucune information pertinente sur le mouvement ondulatoire de l’équation (21). La solution u(x, t) = 0 signifie
qu’il n’existe aucun mouvement ondulatoire ! c’est une solution dite triviale. Ce qui nous amène à dire :
Cours complet est disponible sur mon site web : [Link]

Toutes les solutions physiquement acceptables sont solutions de l’équation (21), mais toutes les solutions de
(21) ne sont pas physiquement acceptables.

o Cas où k > 0, posons k = β 2 , β ∈ R

00 0
Posons ψ = λ2 (donc ψ = λ1 et ψ = λ0 = 1) et écrivons le polynôme caractéristique de l’équation (25)
qui devient :

λ2 − β 2 = 0 ⇒ λ1,2 = ±β (28)
La solution générale de l’équation (25), pour k positif, prend la forme :

ψ(x) = c1 eλ1 x + c2 eλ2 x = c1 eβ x + c2 e−β x (29)

Appliquons les conditions aux limites de (22), il vient :

ψ(x = 0) = c1 + c2 = 0 ⇒ c2 = −c1

ψ(x = l) = c1 eβ l + c2 e−β l = c1 (eβ l + e−β l ) = 0 ⇒ c1 = 0 ⇒ c2 = 0

De façon analogue que précédemment, nous obtenons une solution triviale !. Intéressons-nous désormais au
dernier cas.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 62

 M. Samir KENOUCHE

o Cas où k < 0, posons k = −β 2 , β ∈ R

MASTER PHYSIQUE - Toutes options confondues

Écrivons le polynôme caractéristique de l’équation (25) qui devient :

λ2 + β 2 = 0 ⇒ λ21,2 = −β 2 ⇒ λ21,2 = j 2 β 2 ⇒ λ1,2 = ±j

p
β (30)
La solution générale de l’équation (25), pour k négatif, prend la forme :

ψ(x) = c1 eλ1 x + c2 eλ2 x = c1 ej β x + c2 e−j β x (31)

Afin de simplifier les calculs, écrivons l’équation (31) sous forme d’une combinaison de fonctions sinusoı̈dales.
Utilisons pour cela la formule d’Euler :

e±j θ = cos(θ) ± j sin(θ) (32)

A partir de l’équation (31) :

ψ(x) = c1 cos(β x) + c1 j sin(β x) + c2 cos(β x) − c2 j sin(β x)

= (c1 + c2 ) cos(β x) + (c1 j − j c2 ) sin(β x)
| {z } | {z }
cα ∈R cβ ∈C

⇒ ψ(x) = cα cos(β x) + cβ sin(β x) (33)

Comme précédemment, les constantes cα et cβ sont déterminées à partir des conditions aux limites (22), soit :

ψ(x = 0) = cα = 0
ψ(x = l) = 0 = cα cos(β l) + cβ sin(β l)
Cours complet est disponible sur mon site web : [Link]

⇒ ψ(x = l) = cβ sin(β l) = 0 (34)

Cette équation est nulle dans deux cas de figure. D’abord si cβ = 0 dans ce cas il en résulte cα = cβ =
0 ⇒ ψ(x) = 0 c’est une solution triviale qui n’est pas intéressante d’un point de vue physique. Ensuite, la
deuxième condition si :

cβ 6= 0 ⇒ sin(β l) = 0 ⇒ β l = n π, avec n ∈ N∗ (35)

Ainsi, la solution de l’équation (25) pour k < 0 prend la forme :

nπx
 
ψ(x) = cβ sin avec n ∈ N∗ (36)
l
Cette solution décrit l’amplitude spatiale de l’onde de matière en fonction de la position. Résolvons désormais
00
f (t) − k v 2 f (t) = 0 pour k = −β 2 , avec β ∈ R soit :
00
f (t) + β 2 v 2 f (t) = 0 (37)
De manière analogue que précédemment, au moyen du polynôme caractéristique nous obtenons la solution
générale :

f (t) = c1 eλ1 t + c2 eλ2 t = c1 ej β v t + c2 e−j β v t (38)

Chapitre3 Équations aux dérivées partielles et fonctions de Green 63

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Les termes de la solution (38) sont oscillatoires, par conséquent la quantité v β doit forcément valoir les
dimensions d’une pulsation w soit :

⇒ f (t) = c1 ej ω t + c2 e−j ω t (39)

Qui peut s’écrire également sous la forme équivalente :

f (t) = A cos(ω t + φ) (40)

Tenant compte des solutions (36) et (40), la solution de l’équation (21) devient :
∞
nπx
X  
u(x, t) = An cos(ωn t + φn ) sin (41)
n=1
l
Ce n’est pas cette solution qui nous intéresse dans ce cas précis. Elle est donnée à titre informatif. Le but
est d’obtenir une solution générale de f (t) une fois la nature (positive, négative ou nulle) de la constate de
séparation k est connue. Désormais nous pouvons écrire :

u(x, t) = ψ(x) A cos(ω t + φ) (42)

En substituant (42) dans (23) :

1 d2 ψ(x) −A w2 cos(ω t + φ)
= (43)
ψ(x) dx2 v 2 A cos(ω t + φ)

w2 00
⇒ ψ (x) +
ψ(x) = 0 (44)
v2
Par ailleurs, l’énergie totale d’une particule est la somme des parties cinétique et potentielle soit :
Cours complet est disponible sur mon site web : [Link]

p2
E= + V (x) (45)
2m
En tirant la quantité de mouvement p :

p = {2m[E − V (x)]}1/2 (46)

En utilisant la formule de de Broglie pour la longueur d’onde :
h h
λ= = (47)
p {2m[E − V (x)]}1/2
Le terme ω 2 /v 2 peut être réécrit en fonction de λ, nous rappelons que ω = 2πν et νλ = v.

ω2 4π 2 ν 2 4π 2 2m[E − V (x)]
= = = (48)
v2 v2 λ2 ~2
En substituant ce dernier résultat dans l’équation (44), nous obtenons la fameuse équation de Schrödinger
indépendante du temps :

d2 ψ(x) 2m
+ 2 [E − V (x)]ψ(x) = 0 (49)
dx2 ~
qui est presque toujours écrite sous la forme :

~2 d2 ψ(x)
− + V (x)ψ(x) = Eψ(x) (50)
2m dx2

Chapitre3 Équations aux dérivées partielles et fonctions de Green 64

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Cette équation unidimensionnelle à une seule particule peut facilement être étendue au cas tridimensionnel :

~2 2
∇ ψ(r) + V (r)ψ(r) = Eψ(r)
− (51)
2m
Cette équation peut également traiter un problème à deux corps en remplaçant m par une masse réduite
m2 m1
µ= . Soulignons que Schrödinger a d’abord présenté son équation indépendante du temps, ensuite
m1 + m2
il a postulé l’équation plus générale dépendante du temps.

2) Équation dépendante du temps: examinons désormais l’équation de Schrödinger dépendante du

temps. Dans la section précédente, l’équation de Schrödinger indépendante du temps d’une particule a été
déterminée à partir de l’équation des ondes classique et de la relation de de Broglie. En revanche, l’équation de
Schrödinger dépendante du temps ne peut être obtenue au moyen de méthodes élémentaires et est généralement
donnée comme postulat de la mécanique quantique. L’équation de Schrödinger dépendante du temps à une
seule particule est la suivante :

∂ψ(r, t) ~2 2
j~ =− ∇ ψ(r, t) + V (r)ψ(r, t) (52)
∂t 2m
Où V est supposé être une fonction réelle et représente l’énergie potentielle à laquelle est soumise la particule.
Notons que l’équation (52) ne tient pas encore compte des effets de spin ou relativistes. Bien entendu, l’équation
dépendante du temps peut être utilisée afin d’établir l’équation indépendante du temps. Si nous écrivons la
fonction d’onde comme un produit de deux fonctions spatiale et temporelle, ψ(r, t) = ψ(r)f (t), alors l’équation
(52) devient :
" #
j~ df 1 ~2 2
= ∇ + V (r) ψ(r) (53)
f (t) dt ψ(r) 2m
Puisque le terme de gauche de l’équation est dépendant uniquement du temps et le terme de droite dépend
uniquement de l’espace, l’égalité de l’équation (53) est satisfaite dans le cas où les deux termes sont égaux
Cours complet est disponible sur mon site web : [Link]

à la même constante. Si nous désignons cette constante E (puisque le côté droit doit clairement avoir les
dimensions de l’énergie), nous en obtenons deux équations différentielles ordinaires, à savoir :
1 df (t) jE
=− (54)
f (t) dt ~
et

~2 2
−
∇ ψ(r) + V (r)ψ(r) = Eψ(r) (55)
2m
Cette dernière équation est celle de Schrödinger indépendante du temps. La solution de (54) est :

f (t) = e−j E t/~ avec Re[e−j E t/~ ] = cos(ω t) (56)

Nous retrouvons le résultat de f (t) écrit pour le cas de l’équation de Schrödinger indépendante du temps.
L’Hamiltonien de l’équation (55) est un opérateur hermitien et les valeurs propres d’un opérateur hermitien
doivent être réelles, donc la constante E est réelle. Cela signifie que les solutions f (t) sont purement oscillatoires,
rappelons la formule d’Euler e±iθ = cos(θ) ± i sin(θ). Par voie de conséquence si :

ψ(r, t) = ψ(r) e−j E t/~ (57)

alors la fonction d’onde totale ψ(r, t) diffère de ψ(r) uniquement par un facteur de phase d’amplitude
constante. Cela a des conséquences intéressantes. Tout d’abord, la quantité |ψ(r, t)|2 est indépendante du
temps, car nous pouvons

|ψ(r, t)|2 = ψ ∗ (r, t) ψ(r, t) = ejEt/~ ψ ∗ (r, t) e−jEt/~ ψ(r, t) = ψ ∗ (r)ψ(r)

Chapitre3 Équations aux dérivées partielles et fonctions de Green 65

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Deuxièmement, la valeur attendue pour tout opérateur indépendant du temps est également indépendante
du temps, si ψ(r, t) satisfait l’équation (57). Par le même raisonnement :
Z Z
∗
hAi = ψ (r, t) Â ψ(r, t) = ψ ∗ (r) Â ψ(r)

Pour ces raisons, les fonctions d’onde de la forme (57) sont appelées états stationnaires. L’équation (57)
représente une solution particulière de l’équation (52). La solution générale de l’équation (52) serait une
combinaison linéaire de ces solutions particulières :

ci ψi (r) e−jEi t/~

X
Ψ(r, t) =
i

B. Solution exacte de l’équation de Schrödinger

La résolution exacte de l’équation de Schrödinger en coordonnées cartésiennes est inextricable pour l’atome
d’hydrogène ou les ions hydrogénoı̈des 2 (He+ , Li2+ , · · · etc) à cause de la non séparabilité des variables. Cette
difficulté est levée si l’on considère un système de coordonnées sphériques dont les variables sont séparables. Les
coordonnées sphériques facilitent grandement la résolution exacte de l’équation de Schrödinger pour l’atome
d’hydrogène et les ions hydrogénoı̈des. Ainsi avant de rentrer dans le vif du sujet, nous commencerons par
résumer les principaux résultats d’un tel système de coordonnées.

z
dr

r sinθ
dϕ
e'
e
rdθ

θ r
Cours complet est disponible sur mon site web : [Link]

dθ

ϕ y

dϕ
r sinθ dϕ
x

Figure 1: Représentation des coordonnées sphériques

Il faut bien garder à l’esprit que lors de l’intégration volumique, il est judicieux de savoir que l’élément de
volume en question dépend du système de coordonnées. L’élément de volume forme un parallélépipède dont les
arêtes quantifient les déplacements élémentaires obtenus lorsque l’on fait varier une seule des trois coordonnées.
Dans le système de coordonnées sphériques, les déplacements élémentaires s’écrivent selon dv = du dv dw. Le
déplacement de l’électron dans la direction radiale (r) conduit à une variation élémentaire du = dr. En gardant
la distance radiale et l’angle azimutal (φ) fixes, l’abscisse curviligne générée par la variation de l’angle θ devient
v = r θ => dv = r dθ. Avec un raisonnement analogue, nous obtenons pour le dernier déplacement élémentaire
dw = r sin(θ) dφ. Ce qui donne comme élément de volume :

dv = r2 sin(θ) dr dθ dφ (58)

2. Atomes ayant une structure électronique semblable à celle de l’atome d’hydrogène.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 66

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Il convient de remarquer également que cet élément de volume n’est pas constant, car il dépend de la distance
radiale et de l’angle θ. L’écriture de l’équation de Schrödinger d’un tel système donne :
!
~2 2
− ∇ +V ψ =Eψ (59)
2m
Où ∇2 est l’opérateur Laplacien et V est l’énergie potentielle électrostatique ou Coulombienne qui est donnée
par :

1 Z e2
V (r) = − (60)
4 π 0 r
Où 0 est la permittivité du vide (pas besoin d’une permittivité relative car l’espace à l’intérieur de l’atome
est ”vide”), les charges e et Ze sont respectivement celles de l’électron et du noyau, pour l’hydrogène et
les ions hydrogénoı̈des le nombre d’électron Z = 1. La distance radiale, r, décrit l’éloignement de l’électron
par rapport au noyau. L’énergie potentielle Coulombienne est inversement proportionnelle à la distance entre
l’électron et le noyau et ne dépend d’aucun angle. Un tel potentiel est appelé potentiel central. La solution
exacte de l’équation de Schrödinger pour l’atome d’hydrogène et les ions hydrogénoı̈des est obtenue sous la
forme :

ψ(r, θ, φ) = Rn,l (r) × Ylm (θ, φ) (61)

| {z } | {z }
taille de l’orbitale forme de l’orbitale

En utilisant le Laplacien en coordonnées sphériques, l’équation de Schrödinger devient :

" #
~2 1 ∂ ∂ 1 ∂ ∂ 1 ∂2
   
− 2
r2 + 2 sin θ + 2 2 ψ + V (r)ψ = E ψ (62)
2m r ∂r ∂r r sin θ ∂θ ∂θ r sin θ ∂φ2
| {z }
∇2

1 ∂ ∂ψ 1 ∂ ∂ψ 1 ∂2ψ 2 m
   
Cours complet est disponible sur mon site web : [Link]

⇒ 2 r2 + 2 sin θ + 2 + 2
r ∂r ∂r r sin θ ∂θ ∂θ r2 sin θ ∂φ2 ~
!
Z e2
× E+ ψ=0 (63)
4π 0 r
En utilisant la méthode de séparation des variables, nous considérons une solution ψn,l,m (r, θ, φ) s’écrivant
comme un produit d’une fonction radiale Rn,l (r) et d’une fonction angulaire Yl,m (θ, φ) :

ψn,l,m (r, θ, φ) = Rn,l (r) × Yl,m (θ, φ) (64)

Ce qui donne :

Y ∂R 2 ∂R R ∂ ∂Y R ∂2Y 2m
   
⇒ 2 r + 2 sin θ + 2 2 2
+ 2
r ∂r ∂r r sin θ ∂θ ∂θ r sin θ ∂φ ~
!
Z e2
× E+ RY = 0 (65)
4π 0 r
Désormais nous multiplions par r2 et divisons par R Y afin de séparer la variable radiale et les variables
angulaires :

1 d dR 1 ∂ ∂Y 1 ∂2Y 2 m r2
   
⇒ r2 + sin θ + +
R dr dr Y sin θ ∂θ ∂θ Y sin2 θ ∂φ2 ~2
!
Z e2
× E+ =0 (66)
4π 0 r

Chapitre3 Équations aux dérivées partielles et fonctions de Green 67

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

!
1 d dR 2 m r2 Z e2
 
⇒ r2 + E+ +
R dr dr ~2 4π 0 r
| {z }
Partie radiale

1 ∂ ∂Y 1 ∂2Y
 
sin θ + =0 (67)
Y sin θ ∂θ ∂θ Y sin2 θ ∂φ2
| {z }
Partie angulaire

Les deux parties s’annulent dans le cas où les deux termes (radial et angulaire) sont égaux à la même
constante mais de singe opposé. La constante choisie est connue sous le nom de constante de séparation,
notons cette constante K. Ainsi nous obtenons les deux équations différentielles suivantes :
!
d dR 2 m r2 Z e2
 
⇒ r2 + E+ R−KR=0 (68)
dr dr ~2 4π 0 r
| {z }
Partie radiale

1 ∂ ∂Y 1 ∂2Y
 
⇒ sin θ + +KY =0 (69)
sin θ ∂θ ∂θ sin2 θ ∂φ2
| {z }
Partie angulaire

1) Résolution de la partie angulaire: La partie angulaire contient encore des termes dépendant à la
fois de θ et φ. Une autre séparation des variables est nécessaire. Remplaçons la fonction angulaire Y (θ, φ) par
le produit :

Y (θ, φ) = f (θ) × g(φ) (70)

g d df f d2 g
 
⇒ sin θ + +Kfg =0 (71)
Cours complet est disponible sur mon site web : [Link]

sin θ dθ dθ sin2 θ dφ2

En isolant les deux variables :

sin θ d df 1 d2 g
 
⇒ sin θ + K sin2 θ + =0 (72)
f dθ dθ g dφ2
De la même façon que précédemment notons B la constante de séparation, nous obtenons les deux équations
différentielles suivantes :
sin θ d df
 
⇒ sin θ + K sin2 θ −B = 0 (73)
f dθ dθ
| {z }
Partie polaire

1 d2 g
⇒ +B = 0 (74)
g dφ2
| {z }
Partie azimutale

La solution générale de la partie azimutale est donnée par :

g(φ) = c1 eλ1 φ + c2 eλ2 φ (75)

La solution de la partie azimutale s’obtient en écrivant le polynôme caractéristique de l’équation différentielle
en question :
√
λ2 + B = 0 ⇔ λ21,2 = j 2 B ⇒ λ1,2 = ± j B (76)

Chapitre3 Équations aux dérivées partielles et fonctions de Green 68

 M. Samir KENOUCHE

Il est clair que la constante B doit être positive. Notons cette constante m2 (donc B = m2 ) ⇒ λ1 = +j m

MASTER PHYSIQUE - Toutes options confondues

et λ2 = −j m. Ainsi la solution générale prend la forme :

g(φ) = c1 ej m φ + c2 e−j m φ (77)

L’angle φ est l’azimut, c’est-à-dire que si nous considérons l’atome comme un globe, alors c’est la longi-
tude de la position de l’électron. Nous pouvons choisir le méridien de Greenwich de l’atome d’une manière
mathématiquement commode en fixant c2 = 0. En terminologie quantique, m est appelé un nombre quantique 3
car il limite les valeurs possibles de la fonction d’onde (et donc des observables) à des multiples entiers.

gm (φ) = c1 ej m φ (78)
L’indice m est ajouté à gm (φ) car il est désormais clair qu’il existe autant de solutions qu’il existe des valeurs
autorisées de m. La condition de périodicité de l’angle φ impose :

gm (φ) = c1 ej m φ = c1 ej m (φ+2π) = c1 ej m φ ej m 2π ⇒ ej m 2π = 1 (79)

⇒ m = 0, ±1, ±2, ±3, · · · (80)

Considérons désormais la partie polaire dont l’équation différentielle s’écrit :
sin θ d df
 
sin θ + K sin2 θ − m2 = 0 (81)
f dθ dθ
Qui se réarrange :
!
1 d df m2
 
sin θ + K − f =0 (82)
sin θ dθ dθ sin2 θ
d dx d d
Cours complet est disponible sur mon site web : [Link]

Posons x = cos θ ⇒ = = − sin θ alors l’équation précédente devient :

dθ dθ dx dx
!
1 d df m2
 
(− sin θ) sin θ(− sin θ) + K − f =0 (83)
sin θ dx dx sin2 θ
En exploitant la relation trigonométrique sin2 θ + cos2 θ = 1 ⇒ sin2 θ = 1 − cos2 θ = 1 − x2 . Il en résulte :
!
d df m2
 
(1 − x2 ) + K − f =0 (84)
dx dx 1 − x2
En différentiant le premier terme de l’équation, nous obtenons l’expression finale :
!
d2 f
2 df m2
(1 − x ) 2 − 2 x + K − f =0 (85)
dx dx 1 − x2
La formule finale de cette équation différentielle (n’oublions pas que nous avons déjà posé x = cos θ) est
bien connue dans la littérature mathématique. Les solutions de cette équation sont obtenues pour K = l(l + 1)
et sont connues sous le nom de polynômes de Legendre :
s
m (2 l + 1)! (l − m)!
fl,m (θ) = (−1) × Pl,m (cos θ) pour m ≥ 0 et n ≥ l + 1 (86)
4 π (l + m)!
Avec,
q dm
Pl,m (cos θ) = (−1)m (1 − cos2 θ)m Pl (cos θ) (87)
d cosm θ
3. C’est le nombre quantique azimutal

Chapitre3 Équations aux dérivées partielles et fonctions de Green 69

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Où,

(−1)l dl
Pl (cos θ) = (1 − cos2 θ)l (88)
2l ! d cosl θ
Les polynômes associés de Legendre Pl,m (cos θ) sont construits à partir des polynômes de Legendre Pl (cos θ).
La solution finale 4 de la partie angulaire s’écrit comme le produit des solutions des parties polaire et azimutale
soit :

Yl,m (θ, φ) = fl,m (θ) × ej m φ (89)

Dans la littérature mathématique les solutions Yl,m (θ, φ) sont appelées harmoniques sphériques 5 . Ces fonc-
tions sont particulièrement utiles pour résoudre des problèmes invariants par rotation. Pour résumer cette
section, en utilisant la méthode de séparation des variables grâce aux coordonnées sphériques, nous avons
séparé l’équation angulaire en deux parties azimutale (dépendant de φ) et polaire (dépendant de θ). Les
solutions de l’équation de l’angle azimutal sont des exponentielles incluant le nombre quantique magnétique m
comme argument. Les solutions de l’équation de l’angle polaire sont les polynômes associés de Legendre, qui
sont différents pour chaque choix du nombre quantique azimutal l et de nombre quantique magnétique m. Les
deux nombres quantiques sont introduits dans les équations différentielles respectives en tant que constantes
de séparation.

2) Résolution de la partie radiale: dans ce qui suit, nous développerons une approche étape par étape
afin de résoudre la partie radiale de l’équation de Schrödinger pour l’atome d’hydrogène et les hydrogénoı̈des.
Les énergies propres négatives de l’Hamiltonien sont recherchées comme solution, car elles représentent les
états liés de l’atome. Nous avons déjà obtenu pour la partie radiale l’expression suivante :

d dR 2 m r2
 
r2 − (V (r) − E) R = l(l + 1) R (90)
dr dr ~2
Où la constante de séparation K = l(l + 1) a été obtenue pour la partie angulaire. Ci-dessous les étapes
Cours complet est disponible sur mon site web : [Link]

détaillées de la résolution de la partie radiale.

¶ + : Nous devons d’abord simplifier l’équation radiale pour faciliter la résolution de l’équation différentielle.
Les sous-étapes suivantes utilisent la technique des substitutions pour créer une équation différentielle résoluble.
u(r) dR(r) du(r) 1
 
Posons u(r) = r R(r) ⇒ R(r) = ⇒ = r −u × 2
r dr dr r
" #
2 2 2 2
−~ d u(r) e ~ l(l + 1)
+ − + u=Eu (91)
2 m dr2 4π 0 r 2 m r2
−2 m E
Posons γ 2 = il vient :
~
" #
1 d2 u(r) m e2 1 l(l + 1)
= 1− × + u (92)
γ 2 dr2 2π 0 ~2 γ γ r γ 2 r2
m e2
Posons ρ = γ r et ρ0 = nous obtenons :
2π 0 ~2 γ
d2 u(r) ρ0 l(l + 1)
 
= 1 − + u (93)
dρ2 ρ ρ2
· + : Maintenant que l’équation est sous une forme appropriée pour la solution. Cette étape consiste à
identifier les points singuliers. Il existe des points singuliers où la fonction d’onde doit tendre vers zéro. Dans
∗
4. Pour m < 0 nous avons Yl,−m (θ, φ) = Yl,m (θ, φ)
5. Ces harmoniques sphériques sont normalisées, ce qui explique la disparition de la constante d’intégration c1 obtenue pour
la partie azimutale

Chapitre3 Équations aux dérivées partielles et fonctions de Green 70

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

ce cas, la fonction d’onde doit s’annuler au centre de l’atome, donc pour r = 0. Elle doit également s’annuler à
une distance relativement ”grande” du noyau, prise comme r −→ ∞. Chaque point singulier doit être considéré
individuellement. A ce stade nous étudierons donc le comportement asymptotique de u(ρ) :
ρ0 l(l + 1)
– Pour ρ −→ ∞ ⇒ −→ 0 et −→ 0. Comme la distance r de l’atome va à l’infini, ces
ρ ρ2
termes tendent vers zéro et sont donc sans importance pour cette partie de la solution. Ainsi, l’équation
différentielle à résoudre, sous la condition d’une distance infinie, devient :

d2 u(r)
⇒ =u (94)
dρ2
La solution générale de l’équation différentielle ainsi obtenue est : u(ρ) = A1 e−ρ + A2 eρ . Le deuxième
terme de la solution eρ est refusé d’un point de vue physique. Car −→ ∞ ⇒ eρ −→ ∞ or l’électron
possède une limite spatiale par rapport au noyau. Il en ressort :

=⇒ u(ρ) ' e−ρ (95)

l(l + 1) ρ0
– Le deuxième point singulier c’est quand r −→ 0 (donc au centre de l’atome) ⇒ 2
>> et
ρ ρ
l(l + 1)
>> 1. Par conséquent, à ce deuxième point singulier, où r tend vers zéro, l’équation différentielle
ρ2
devient :

d2 u(r) l(l + 1)
2
⇒ = u (96)
dρ ρ2
La solution générale de cette équation différentielle est de la forme :

=⇒ u(ρ) ' ρl+1 (97)

En combinant les solutions obtenues pour les deux points singuliers (95) et (97), la solution complète de
Cours complet est disponible sur mon site web : [Link]

l’équation (93) prend la forme :

=⇒ u(ρ) = ρl+1 e−ρ y(ρ) (98)

Où y(ρ) est un polynôme exprimé sous forme :
∞
X
y(ρ) = cj ρj (99)
j=0

Nous nous sommes servi du comportement asymptotique de u(ρ) pour trouver l’expression de u(ρ) pour
0 < ρ < ∞.

¸ + : Après avoir obtenu une forme générale de la solution complète. Nous devons maintenant trouver une
équation pour la partie polynomiale y(ρ) de cette solution complète. Pour cela, calculons les dérivées première
et second de l’équation (98) :
du(ρ) dy(ρ)
 
= ρl e−ρ (ρ + 1 − 1) y + ρ (100)
dρ dρ
et,

du(ρ) l(l + 1) dy(ρ) d2 y(ρ)

 
= ρl e−ρ (−2ρ − 2 + l + ) y + 2 (l + 1 − ρ) +ρ (101)
dρ l dρ dρ2
En substituant (100) et (101) dans (93) nous obtenons :

d2 y(ρ) dy(ρ)
ρ 2
+ 2 (l + 1 − ρ) + [ρ0 − 2 (l + 1)] y = 0 (102)
dρ dρ
Chapitre3 Équations aux dérivées partielles et fonctions de Green 71
 M. Samir KENOUCHE

En utilisant les séries entières 6 , nous chercherons des solutions de y(ρ) de la forme :

MASTER PHYSIQUE - Toutes options confondues

∞
X
y(ρ) = cj ρj (103)
j=0

Où les inconnus sont les coefficients cj . Les dérivées de (103) se calculent selon :
∞ ∞
dy(ρ) X X
= j cj ρj−1 = j + 1 cj+1 ρj (104)
dρ j=0 j=0

et,
∞
d2 y(ρ) X
= j (j + 1) cj+1 ρj−1 (105)
dρ2 j=0

Substituons (104) et (105) dans l’équation (102) :

∞
X ∞
X ∞
X
ρ j (j + 1) cj+1 ρj−1 + 2(l + 1) (j + 1) cj+1 ρj − 2ρ j cj ρj−1 , +
j=0 j=0 j=0
∞
X
[ρ0 − 2(l + 1)] cj ρj = 0 (106)
j=0

Qui se simplifie :
∞
X ∞
X ∞
X
j (j + 1) cj+1 ρj + 2(l + 1) (j + 1) cj+1 ρj − 2 j cj ρj + [ρ0 − 2(l + 1)]
j=0 j=0 j=0
∞
X
× cj ρj = 0 (107)
Cours complet est disponible sur mon site web : [Link]

j=0

Nous n’avons pas encore fini avec le polynôme en question, car nous devons déterminer la relation de
récurrence pour ses coefficients. Nous devons déterminer cette relation non seulement pour savoir comment le
polynôme sera généré, mais aussi pour déterminer les limites de la sommation de la série. Le polynôme (107)
vaut zéro si et seulement si cj = 0, soit :

j (j + 1) cj+1 + 2(l + 1) (j + 1) cj+1 − 2 j cj + [ρ0 − 2(l + 1)] cj = 0 (108)

Nous obtenons la formule de récurrence suivante :
2 (j + l + 1) − ρ0 2 (j + l + 1) − ρ0
cj+1 = cj = cj (109)
j (j + 1) + 2 (l + 1) (j + 1) (j + 1) [j + 2 (l + 1)]
Cette formule de récurrence décrit le comportement des coefficients de la série (du polynôme). Le polynôme
y(ρ) doit tendre vers zéro, donc le comportement du polynôme doit être examiné lorsque j −→ ∞. Cela fera
l’objet de la prochaine étape.

¹ + : Dans cette dernière étape de résolution de l’équation radiale, nous examinons le polynôme pour
déterminer s’il est fini. Sinon, nous déterminons quelle condition est nécessaire pour le rendre fini. Cette
condition de finitude produit un nombre quantique qui caractérise l’état du système et sert à quantifier

6. Le lecteur est renvoyé aux références :

Griffiths, Introduction to Quantum Mechanics, Prentice Hall, Englewood Cliffs, New Jersey, (1995), pp. 134-141. ET Cohen-
Tannoudji, Diu, and Laloe, Quantum Mechanics, John Wiley & Sons, New York, (1977), pp. 794-797.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 72

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

les énergies d’état lié de l’atome. En théorie, la formule (109) peut se développer à l’infinie, étudions son
comportement quand j −→ ∞ (c’est-à-dire pour j >> 1) :
2j 2
j >> 1 ⇒ cj+1 ' cj ⇒ cj+1 ' cj (110)
j (j + 1) (j + 1)
22 2j
⇒ cj+2 ' cj , ······ , cj ' c0 (111)
(j + 1) (j + 2) j!
En substituant (111) dans (103) nous obtenons :
∞
X 2j
y(ρ) = c0 ρj (112)
j=0
j!
| {z }
Devlop. Taylor de e2ρ

En substituant (112) dans (98) nous obtenons :

=⇒ u(ρ) = c0 ρl+1 e−ρ e2ρ ⇒ lim u(ρ) −→ ∞ (113)

ρ→∞

La solution (113) diverge lorsque r est très grand. Ce comportement de la solution n’est pas acceptable car
l’électron possède une limite spatiale finie par rapport à sa distance du noyau. Par conséquent, la série doit
être tronquée à un nombre particulier, afin de forcer le polynôme d’avoir un comportement correct. Afin de
déterminer le rang où doit se produire la troncature, nous fixons le coefficient du polynôme y(ρ) égal à zéro
à un nombre maximum, forçant la terminaison de la solution à de grandes distances du noyau. A partir de la
formule de récurrence (109), nous avons :

cjmax +1 = 0 ⇒ 2 (jmax + l + 1) −ρ0 = 0 ⇒ ρ0 = 2 n (114)

| {z }
n
Où n est le nombre quantique principal. La nouvelle formule de récurrence s’obtient :
Cours complet est disponible sur mon site web : [Link]

2 (j + l + 1) − 2 n
cj+1 = cj (115)
(j + 1) [j + 2 (l + 1)]
1
Appliquons cette relation par exemple pour n = 3 et l = 1 ⇒ jmax = 1 ⇒ c1 = − c0 :
2
l 1 1 1
⇒ yn (ρ) = y3 (ρ) = c0 − c0 ρ = c0 (1 − ρ) (116)
2 2
On peut continuer ce processus à l’infini, mais on cherche une solution analytique à cette équation. La forme
asymptotiquement suggérée nous donne un point de départ pour chercher la solution finale. Tenant compte
de cette forme, nous soupçonnons une solution de la forme :

w(ρ) = ρl+1 × e−ρ (117)

0
w = −ρl+1 e−ρ + (l + 1) ρl e−ρ (118)
00 −ρ l+1 l −ρ l−1 l −ρ
w =e ρ − (l + 1) ρ e + l (l + 1) ρ − (l + 1) ρ e (119)
Pour l = 1, nous obtenons :

00
w = e−ρ ρ2 − 2 ρ e−ρ + 2 − 2 ρ e−ρ (120)
00 −ρ 2
w =e [ρ − 2 ρ + 2 − 2 ρ] (121)
00 −ρ 2
w =e [ρ − 2 ρ + 2 − 2 ρ] (122)
00
eρ w = [ρ2 − 4 ρ + 2] (123)

Chapitre3 Équations aux dérivées partielles et fonctions de Green 73

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

En dérivant la dernière équation nous obtenons :
00
d[eρ w ] 1
= 2ρ − 4 = ρ − 1 (124)
dρ 2
En comparant (116) à (124) nous déduisons :
00
d[eρ w (ρ)]
y31 (ρ) = (−1)3 (125)
dρ
Cette dernière relation s’écrit également sous la forme :
3 " 2 #
d d
 
ρ −ρ l+1
y31 (ρ) = (−1) 3
e e ρ (126)
dρ dρ
Le deuxième terme de l’équation (126) n’est autre que le polynômes associès de Laguerre d’ordre n = 3,
notée L23 (2 ρ). En définitif, la généralisation de ce résultat est immédiate et nous obtenons :

y(ρ) = L2l+1
n−l−1 (2 ρ) ∀ n ≥ l + 1 (127)
Nous rappelons que les polynômes associés de Laguerre sont définis par les relations suivantes :
q
dp ρ d

Lqp p
e−ρ ρq
 
= (−1) e (128)
dρp dρ
En combinant (98) et (127) nous obtenons la solution finale 7 :

u(ρ) = (2 ρ)l+1 e−ρ L2l+1

n−l−1 (2 ρ) ∀ n ≥ l + 1 (129)
−2 m E
Au début de la résolution de la partie radiale, nous avons déjà posé : γ 2 = , ρ = γ r et ρ0 =
~
m e2
Cours complet est disponible sur mon site web : [Link]

= 2n.
2 π 0 ~2 γ
~2 γ 2 m e4 m e4
⇒ E= =− 2 2 2 2 =− 2 2 2 (130)
2m 8 π 0 ~ γ0 8 π 0 ~ (2 n)2
" #2
m e2 1
⇒ En = − 2 × (131)
2~ 4 π 0 n2
D’un autre côté :

−m E 2 m2 e4 m e2 1 1
⇒ γ2 = = ⇒ γ= × ⇒ γ=
~2 8 π 2 20 ~2 (2n)2 4 π 0 ~ n a0 n
| {z }
a0
−10
Avec a0 = 0.53 10 m est le rayon de Bohr. Revenons maintenant à la fonction radiale initiale R(r) par le
changement de variable que nous avons réalisé au début de notre résolution soit :
u(r) r
u(r) = r R(r) ⇒ R(r) = Où ρ = γ r = (132)
r a0 n
r
 
l+1 −
1 2r 2r
  
Rn,l (r) = e a0 n L2l+1 (133)
n−l−1
r a0 n a0 n

7. Les polynômes associés de Laguerre, produits par la troncature de la série, sont identifiés par deux indices ou nombres
quantiques, n et l. Les solutions physiquement acceptables exigent que n soit supérieur ou égal à l + 1.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 74

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

r
 
l −
1 2r 2r 2r
    
Rn,l (r) = e a0 n L2l+1 (134)
n−l−1
r n a0 a0 n a0 n
r
 
l −
2 2r 2r
   
⇒ Rn,l (r) = e a0 n L2l+1 (135)
n−l−1
n a0 a0 n a0 n
L’équation (135) est la solution finale non normalisée de l’équation radiale. Pour tenir compte de la norma-
lisation :
r
 
l −
2r 2r
  
⇒ Rn,l (r) = Nn,l e a0 n L2l+1 (136)
n−l−1
a0 n a0 n
2
 
Où le facteur est adossé à la constante de normalisation Nn,l qui s’obtient en calculant l’intégrale :
n a0
Z ∞
2
Nn,l Rn,l (r)∗ Rn,l (r) r2 dr = 1 (137)
0

Où r2 provient de l’élément de volume exprimé en coordonnées sphériques. Le calcul de cette constante
étant très laborieux, nous donnons sa valeur :
" 3 #1/2
2 (n − l − 1)!
Nn,l = (138)
n a0 2n [(n + 1)!]3
La solution normalisée de la partie radiale s’écrit alors :
r
 
" #1/2  l
2
3
(n − l − 1)! 2r − 
2r

⇒ Rn,l (r) = e a0 n L2l+1 (139)
n−l−1
n a0 2n [(n + 1)!]3 a0 n a0 n
Cours complet est disponible sur mon site web : [Link]

La solution exacte (valeurs et fonctions propres) de l’équation de schrodinger pour l’atome d’hydrogène (et
les ions hydrogénoı̈des He+ , Li2+ , · · · etc) s’obtient en multipliant les solutions des parties angulaires (89) et
radiale (139) :
r
 
" #1/2  l
2
3
(n − l − 1)! 2r − 
2r

ψn,l,m (r, θ, φ) = e a0 n L2l+1
n−l−1
n a0 2n [(n + 1)!]3 a0 n a0 n
s
m (2 l + 1)! (l − m)!
× (−1) × Pl,m (cos θ) × ej m φ (140)
4 π (l + m)!
| {z }
Ylm (θ,φ)

Ou simplement en occultant la partie angulaire :

r
 
" #1/2  l
2
3
(n − l − 1)! 2r −
⇒ ψn,l,m (r, θ, φ) = e a0 n
n a0 2n [(n + 1)!]3 a0 n
2r
 
× L2l+1
× Ylm (θ, φ)
n−l−1 (141)
a0 n
Le terme exponentiel décroissant supplante le terme polynomial croissant de sorte que la fonction d’onde
globale ψn,l,m (r, θ, φ) tend vers zéro pour les grandes valeurs de r (loin du noyau), c’est ce qui est attendu.
Les valeurs propres sont obtenues avec l’équation (131) :
" #2
m e2 1
⇒ En = − 2 × (142)
2~ 4 π 0 n2

Chapitre3 Équations aux dérivées partielles et fonctions de Green 75

 M. Samir KENOUCHE

Les harmoniques sphériques Ylm (θ, φ), fournissent des informations sur la position de l’électron autour du

MASTER PHYSIQUE - Toutes options confondues

noyau, et la fonction radiale Rn,l (r) décrit l’éloignement de l’électron par rapport au noyau. Comme on peut le
constater, l’équation de Schrödinger requiert trois nombres quantiques (n, l, m) afin de spécifier une fonction
d’onde pour l’électron. Les nombres quantiques fournissent des informations sur la distribution spatiale d’un
électron. Bien que n puisse prendre n’importe quel nombre entier positif non nul, seules certaines valeurs de
l et de m sont autorisées pour une valeur donnée de n. Le nombre quantique principal n indique l’énergie
de l’électron et la distance moyenne d’un électron par rapport au noyau. Plus un électron est proche du
noyau, chargé positivement, plus l’électron est fortement attiré par le noyau comparativement à un électron
plus éloigné dans l’espace. Cela signifie que les électrons ayant une valeur de n plus élevée sont plus faciles à
éliminer d’un atome.

Le deuxième nombre quantique l est appelé nombre quantique azimutal. Ce dernier décrit la forme de la
région de l’espace occupée par un électron, donc la sous-couche considérée. Les valeurs de ce nombre quantique
sont données par n ≥ l + 1. Le troisième nombre quantique, est le nombre quantique magnétique m. Ce
nombre quantique décrit l’orientation de la région dans l’espace occupé par un électron par rapport à un
champ magnétique appliqué. Les valeurs autorisées de m dépendent de la valeur de l selon −l ≤ m ≤ +l.
Chaque combinaison autorisée des trois nombres quantiques fournit une distribution spatiale particulière à
l’électron.

Il est intéressant de comparer les résultats obtenus en résolvant l’équation de Schrödinger avec le modèle de
l’atome d’hydrogène de Bohr. Les valeurs propres (spectre énergétique) sont quasiment identiques. Toutefois,
les modèles de Schrödinger et de Bohr sont différents à bien des égards, notamment en ce qui concerne les deux
points énumérés ci-dessous :

1) Le modèle de Schrödinger n’associe pas d’orbites bien définies pour l’électron. Les fonctions d’onde donnent
seulement la probabilité de trouver l’électron dans l’élément de volume dv à différentes directions (θ et φ)
et distances du noyau (r).
Cours complet est disponible sur mon site web : [Link]

2) Les nombres quantiques apparaissent spontanément lors de la résolution de l’équation de Schrödinger

alors que Bohr a dû postuler l’existence d’états énergétiques quantifiés. Bien que plus complexe, le modèle
de Schrödinger conduit à une meilleure correspondance entre la théorie et l’expérience.

· · · Ouuf c’est terminé · · · dire que l’hydrogène est l’atome le plus simple !

C. Résolution numérique de quelques EDPs

Mis à part certaines EDPs particulières, la grande majorité des EDPs issues de la physique et de la chimie
n’admette pas de solution explicite ou analytique. Il est donc impératif de recourir à la résolution numérique sur
ordinateur pour évaluer qualitativement et quantitativement les solutions. Le principe de base de ces méthodes
de résolution numérique des EDPs, consiste à chercher des valeurs numériques discrètes approchant au mieux
la solution exacte. Le concept le plus important dans cette résolution est celui de discrétisation, marquant la
passage du continu au discret. Dans cette section, on se propose de résoudre numériquement quelques EDPs
régissant des phénomènes bien connus de la physique. La résolution numérique sera conduite en considérant
la méthode des différences finies.

Commençons par chercher les approximations des dérivées première et seconde. Soit f : R2 −→ R une
fonction de classe C 2 (R2 ) et écrivons le développement en séries de Taylor de u(xi + ∆x, tj ) autour du point
(xi , tj ), soit :

∆x ∂u ∆x2 ∂ 2 u ∆xn ∂ n u
u(xi + ∆x, tj ) = u(xi , tj ) + + + ··· + + O(∆xn+1 ) (143)
1! ∂x (xi ,tj ) 2! ∂x2 (xi ,tj ) n! ∂xn (xi ,tj )

Chapitre3 Équations aux dérivées partielles et fonctions de Green 76

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

L’approximation par les différences finies de la dérivée première est obtenue en tronquant la série de Taylor
à l’ordre deux, soit :

∆x ∂u ∆x2 ∂ 2 u
u(xi + ∆x, tj ) = u(xi , tj ) + + + O(∆x3 ) (144)
1! ∂x (xi ,tj ) 2! ∂x2 (xi ,tj )

Après réarrangement, il vient :

∂u u(xi + ∆x, tj ) − u(xi , tj ) ∆x2 ∂ 2 u

= − + O(∆x3 ) (145)
∂x (xi ,tj ) | ∆x
{z } 2! ∂x2 (xi ,tj )
Approx. par DF
| {z }
Erreur de troncature

Cette dernière relation est présentée sous la forme :

∂u u(xi + ∆x, tj ) − u(xi , tj )
= + O(∆x) (146)
∂x (xi ,tj ) | ∆x
{z } | {z }
Erreur de troncature
Approx. par DF

L’écriture de O(∆x) indique que l’approximation de la première dérivée par la formule des différences finie
est d’ordre un par rapport au pas de discrétisation ∆x. Cela signifie concrètement que lorsqu’on divise le
pas de discrétisation ∆x par une constante arbitraire a > 0 implique que l’erreur d’approximation, entre
dérivée exacte et approchée, est divisée par a. Avec un raisonnement analogue, il est possible aussi de calculer
une approximation de la dérivée seconde. Afin d’atteindre cette approximation, nous devons réaliser deux
développement en séries de Taylor (à droite et à gauche), soit :

∆x ∂u ∆x2 ∂ 2 u ∆x3 ∂ 3 u
u(xi + ∆x, tj ) = u(xi , tj ) + + + + O(∆x4 ) (147)
1! ∂x (xi ,tj ) 2! ∂x2 (xi ,tj ) 3! ∂x3 (xi ,tj )

∆x ∂u ∆x2 ∂ 2 u ∆x3 ∂ 3 u
u(xi − ∆x, tj ) = u(xi , tj ) − + − + O(∆x4 ) (148)
1! ∂x (xi ,tj ) 2! ∂x2 (xi ,tj ) 3! ∂x3 (xi ,tj )
Cours complet est disponible sur mon site web : [Link]

En sommant membre par membre les deux dernières relations, nous obtenons l’approximation recherchée :

∂2u u(xi + ∆x, tj ) − 2 u(xi , tj ) + u(xi − ∆x, tj )

= + O(∆x2 ) (149)
∂x2 (xi ,tj ) | ∆x
{z
2
} | {z }
Erreur de troncature
Approx. par DF

Théorème : la solution numérique d’un schéma itératif aux différences finies, d’un problème linéaire aux
valeurs initiales, converge vers la solution exacte si le schéma est consistant et stable.

Définition 1 : une approximation est dite consistante d’ordre p s’il existe une constante arbitraire c > 0
indépendante du pas de discrétisation telle que cette erreur soit majorée par la quantité c ∆xp . Soit u(x, t) ∈
Ω ⊂ R une fonction de classe C 4 (Ω), pour la dérivée d’ordre un, nous avons :

0 u(xi + ∆x, tj ) − u(xi − ∆x, tj ) u(x)(3)

u (x, t) − ≤ max ∆x2 (150)
| 2∆x
{z } x1 ≤x≤x1 +∆x 6
| {z }
Approx. par DF c
| {z }
(u)

Pour la dérivée seconde, il vient :

00 u(xi + ∆x, tj ) − 2 u(xi , tj ) + u(xi − ∆x, tj ) u(x)(4)

u (x, t) − 2
≤ max ∆x2 (151)
| ∆x
{z } x1 ≤x≤x1 +∆x 2
| {z }
Approx. par DF c

Chapitre3 Équations aux dérivées partielles et fonctions de Green 77

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

De manière générale, le schéma numérique des différences finies est dit consistant à l’équation EDP si cette
erreur de troncature tends vers zéro lorsque le pas de discrétisation temporel ∆t et le pas de discrétisation
spatial ∆x tendent indépendamment vers zéro. Autrement dit si,

lim (u) = 0 (152)

∆x≈0,∆t≈0

Définition 1’ : une solution est dite stable 8 si une petite variation des conditions de bord engendre une
faible variation de la solution. En terme mathématique simple cela se traduit par :

∀x ∈ R, ∀t ∈ R+ ⇒ u1 (x, t) − u2 (x, t) ≤ 

Il convient de noter également que l’analyse de la stabilité d’un schéma numérique peut être conduite en
déterminant le facteur d’amplification du schéma itératif.

1) En dimension 1: dans cette section, nous considérons le problème de la corde élastique fixée aux
extrémités x = 0 et x = L telle que L est égale à une unité de longueur. La corde subit des déformations
selon un mode vertical. L’amplitude des déformations est décrite par la fonction u(x), ainsi le problème est
formalisé mathématiquement par :
00




−u (x) = f (x), 0<x<L

P: (153)



u(0) = 0
 u(L) = 0

Autrement dit, l’Eq. (153) est celle régissant la déformation linéaire de la corde élastique. Le second terme
représente la source des déformations. Les conditions aux limites u(0) = 0 et u(L) = 0 traduisent le fait que la
Cours complet est disponible sur mon site web : [Link]

corde ne subit pas de déformations aux extrémités. Il s’agit d’une résolution numérique, donc le calcul d’une
approximation de u(xi ) au point xi . Nous commençons par la discrétisation des abscisses.

h = xi+1 - xi
0 L
x0 x1 xi-1 xi xi+1 xn xn+1

xi = i.h
pour i = 0, 1, ..., n, n+1

Figure 2: Discrétisation du problème aux limites 1 Dim

La discrétisation se fait avec un pas constant h = xi+1 −xi par conséquent tous les points xi sont équidistants.
La solution exacte au point xi , soit u(xi ) est inconnue. Nous cherchons des solutions approchées ui qui sont
des approximations de la solution exacte u(xi ) au point xi .
00
Utilisons une formule des différences finies centrée afin d’approcher u (x), soit :
u(xi−1 ) − 2 u(xi ) + u(xi+1 )
− = f (xi ) + O(h2 ) (154)
h2
8. Cette notion de stabilité consiste à analyser si les perturbations de la solution numérique ne sont pas amplifiées au cours des
itérations. Les calculs sur ordinateur sont déterminés avec une précision finie, sont ainsi sujet à des erreurs d’arrondis. Pendant
un calcul itératif, ces erreurs peuvent être amplifiées.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 78

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

h = xi+1 - xi
0 L
x0 x1 xi-1 xi xi+1 xn xn+1

ui+1

La corde
u(xi+1)

Figure 3: Solution exacte u(xi ) versus solution approchée ui

Le terme O(h2 ) stipule que lorsqu’on divise h par une constante a, l’erreur |u(xi ) − ui | est divisée par a2 .
C’est le résultat du théorème suivant :

Si u(x) est de classe C 4 sur [0 , L], alors ∃ c ∈ R+ telle que ∀ 0 < h < L nous avons :

max |u(xi ) − ui | ≤ c h2
1≤i≤1
1
max |u(xi ) − ui | ≤ max |u(x)(4) |
1≤i≤1 96 0≤x≤L
Ce théorème affirme que l’erreur d’intégration est majorée par la quatrième dérivée de la fonction u(x) et
plus h est petit plus on s’approche de la solution exacte. Écrivons maintenant la même formule des différences
finies pour les approximations ui , il vient :
Cours complet est disponible sur mon site web : [Link]


ui−1 − 2 ui + ui+1

 − = f (xi ), i = 1, 2, ...., n



 h2

u(x0 ) = u0 (155)





 u(L) = un+1


Le terme O(h2 ) n’est pas pris en considération dans les calculs. Le schéma (155) correspond à la résolution
d’un système linéaire.

An (Rn × Rn ) un (Rn ) = fn (Rn ) (156)

Explicitons le schéma (155) pour n = 4 et pour des conditions aux limites u0 = α et un+1 = β. Ces
conditions aux limites signifient qu’à x = 0 la corde subit une déformation constante égale à la valeur α et
à l’autre extrémité x = L, la corde subit aussi une déformation constante égale à la valeur β. Nous aurions
pu prendre par exemple u0 = 0 et un+1 = 0, cela indique que la corde ne subit aucune déformation aux
extrémités. Nous préférons prendre un cas général avec les conditions u0 = α et un+1 = β.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 79

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

−u0 + 2 u1 − u2


 = f (x1 ), i=1
h2

−α + 2 u1 − u2 = h2 f (x1 ), i = 1

 

 

 


 −u1 + 2 u2 − u3 

= f (x2 ), i=2
 
2
 −u1 + 2 u2 − u3 = h f (x2 ), i = 2

 


 h2 
⇒

 −u2 + 2 u3 − u4 
 −u2 + 2 u3 − u4 = h2 f (x3 ), i=3
= f (x3 ), i=3

 

h2

 


 

 
−u3 + 2 u4 − β = h2 f (x4 ), i=4

 
−u3 + 2 u4 − u5





2
= f (x4 ), i=4
h

2 u1 − u2 = h2 f (x1 ) + α, i = 1







2
 −u1 + 2 u2 − u3 = h f (x2 ), i = 2



⇒ (157)
−u2 + 2 u3 − u4 = h2 f (x3 ), i=3









−u3 + 2 u4 = h2 f (x4 ) + β, i=4


En adoptant une notation matricielle, le problème peut s’écrire :

An un = fn (158)
Avec,
 
2 −1 0 0
Cours complet est disponible sur mon site web : [Link]

1 
−1 2 −1 0 

An = 2   ∈ R4 × R4 (159)
h  0 −1 2 −1
0 0 −1 2

 
h2 f (x1 ) + α
 h2 f (x ) 
2
fn =  2  ∈ R4 (160)
 
 h f (x3 ) 
h2 f (x1 ) + β

La matrice An est tridiagonale, symétrique et définie positive. Le vecteur des valeurs de la solution (incon-
nues) aux points xi est donné

 
u1
u 
uh =  2  ∈ R4 (161)
 
u3 
u4
Ce schéma se généralise pour i = {1, 2, ..., n}, selon :

Chapitre3 Équations aux dérivées partielles et fonctions de Green 80

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

 
2 −1 0 ... 0
h2 f (x1 ) + α
   
..  u1
.
−1 . .

−1 2 .   f (x2 )   u2 
     
 .. .. ..  n n
 .
.  ∈ Rn
  .. ∈ Rn

An = 
 0 . . . 0 ∈R ×R fn =  . un =  .
  
 
 .
..  f (xn−1 )  un−1 
    
 .
 . . −1 2 −1

h2 f (xn ) + β un
0 ... 0 −1 2

Exercice ¶ + r

Soit l’équation :
00




−u (x) = f (x) 0≤x≤1

P: (162)



u(0) = 0
 u(1) = 0

Le second terme de l’équation vaut :

2
f (x) = e3 x × (x + 1) (163)

1) Résoudre numériquement pour n = 100 l’équation (162) par la méthode des différences finies.

Voici le script Matlabr

Cours complet est disponible sur mon site web : [Link]

clear␣all␣;␣clc␣;␣close␣all␣;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%␣SAMIR␣KENOUCHE␣-␣RESOLUTION␣NUMERIQUE␣DU␣PROBLEME␣AUX␣LIMITES␣DE␣LA␣CORDE
%␣ELASTIQUE␣:␣-␣u’’(x)␣=␣f(x)␣AVEC␣LES␣CONDITIONS␣u(0)␣=␣0␣et␣u(1)␣=␣0

np␣=␣100␣;␣pas_x␣=␣1/(np+1)␣;␣xi␣=␣0:␣pas_x␣:1␣;␣%␣DISCRETISATION
fx␣=␣@(xi)␣exp(3.*xi.ˆ2).*(xi␣+␣1)␣;

sur_diag␣=␣diag(ones(np␣-␣1,␣1)␣,1)*(-1)␣;␣␣%␣SUR-DIAGONALE
des_diag␣=␣diag(ones(np␣-␣1,␣1)␣,-1)*(-1)␣;␣%␣SOUS-DIAGONALE
in_diag␣=␣diag(ones(np,␣1))*(2)␣;␣␣␣␣␣␣␣␣␣␣␣%␣DIAGONALE␣PRINCIPALE

An␣=␣sur_diag␣+␣des_diag␣+␣in_diag␣;␣␣␣␣␣␣␣␣%␣MATRICE␣An
fn␣=␣fx(xi(2:end-1))␣;␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣%␣SOURCE␣DE␣LA␣DEFORMATION

un␣=␣inv(An)*fn’␣;␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣%␣CALCUL␣DES␣APPROXIMATIONS
un␣=␣[0␣un’␣0]␣;␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣%␣ON␣RAJOUTE␣LES␣CONDITIONS␣AUX␣LIMITES

fig1␣=␣figure(’color’,[1␣1␣1])␣;␣plot(xi,␣un,’o’)␣;
xlabel(’x_i’)␣;␣ylabel(’u_i’)␣;␣title(’SOLUTION␣APPROCHEE’)␣;

Considérons désormais des fonctions c(x) et f (x) continues sur l’intervalle [a, b]. On se propose de résoudre,
par les différences finies, l’équation de la convection-diffusion. Soient α et β deux constantes réelles. Le problème

Chapitre3 Équations aux dérivées partielles et fonctions de Green 81

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

consiste à trouver une fonction u(x) deux fois dérivable sur l’intervalle [a, b] et qui satisfait :

( 00 0
−u (x) + c(x) u (x) = f (x), x ∈ [a, b]
P: (164)
u(a) = 0 et u(b) = 0

Comme précédemment, ce type de problème est dénommé problème aux limites. Cette dénomination provient
du fait que la fonction u(x) doit satisfaire les conditions aux limites, u(a) = 0 et u(b) = 0, posées aux bornes
de l’intervalle [a, b]. Afin de résoudre numériquement (trouver la solution approchée) le système (164), nous
recourrons à la méthode des différences finies. Suivant cette méthode numérique, l’intervalle [a, b] sur lequel nous
cherchons la solution u(x) est discrétisé en n + 1 sous-intervalles équidistants de longueur h avec xi = x0 + i h
et i = 1, 2, 3, ..., n. On cherche alors en chacun de ces points une valeur approchée, notée ui , de u(xi). Ainsi, le
système continu initial est substitué par un système discret. L’idée de base de la méthode des différences finies
consiste à remplacer l’équation différentielle (164) par un système de n équations algébriques. Ce système
d’équations est obtenu en écrivant cette équation différentielle en chaque point de discrétisation xi , et en
00
substituant également à chaque valeur u (x) l’approximation de la dérivée seconde :

00 u(xi−1 ) − 2 u(xi ) + u(xi+1 )

u (x) ≈ + O(h2 ) (165)
h2

0
Et la dérivée u (x) est approchée par

0 u(xi+1 ) − u(xi−1 )
u (x) ≈ + O(h) (166)
2h
Cours complet est disponible sur mon site web : [Link]

Ainsi, l’équation différentielle (164) est réécrite suivant :


u(xi−1 ) − 2 u(xi ) + u(xi+1 ) u(xi+1 ) − u(xi−1 )
− + c(xi ) = f (xi ) i ∈ {1, ..., n}

h 2 2h (167)

u0 = 0 et un = 0

Comme précédemment en adoptant une notation matricielle après quelques réarrangements, le problème
peut s’écrire :

Ah uh = bh (168)
Avec Ah = A1 + A2 ,
 
2 −1 0 ... 0
. .. 
−1 . .

−1 2 . 
 

A1 =  .. .. .. 
(169)
 0 . . . 0

 .
..

 .
 . . −1 2 −1

0 ... 0 −1 2

Chapitre3 Équations aux dérivées partielles et fonctions de Green 82

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

 
0 1 0 ... 0
c(x1 ) 0 ... 0
 
 .. .. 
. .. 
−1 0 1 . .
c(x2 ) . .
  
rh 
 0 .  
× .. .. .. 
A2 = . .. ..   0 . . . 0
 (170)
2 
 .. . . 0   .. ..

 . . −1 0 1

0 ... 0 c(x )
n
0 ... 0 −1 0

 
f (x1 )



 f (x2 )
2

fh = h 
 .. (171)
.


f (xn−1 )
 

f (xn )

Le vecteur des valeurs de la solution (inconnues) aux points xi est donné

 
u1



 u2

uh = 
 ..
 ∈ Rn (172)
  .
un−1 
 

un

On peut mettre en évidence le fait que la matrice Ah est inversible (sous l’hypothèse que la fonction c(x)
Cours complet est disponible sur mon site web : [Link]

≥ 0. La matrice Ah est symétrique définie positive). Les matrices A1 et A2 sont tridiagonales. Afin d’obtenir
la solution discrète, les valeurs du vecteur uh , il va falloir résoudre le système linéaire tridiagonal (168).

Exercice · + r

Soit l’équation de convection-diffusion :

 00 0
 −u (x) + r x u (x) = f (x)

P: (173)
 x ∈ [0, 1]

et u(0) = 0, u(1) = 0

Le second terme de l’équation vaut :

(r2 e(r x) (x − 1))

f (x) = (174)
(1 − er ) + r x

La solution exacte est donnée par :

x − (1 − e(r x) )
funex = (175)
(1 − er )

1) Résoudre numériquement l’équation (173) par la méthode des différences finies.

2) Tracer, sur la même figure, les solutions exacte et numérique pour n=64 et r=1/2.
3) Étudier l’erreur en fonction du nombre de sous-intervalles de discrétisation et du paramètre r.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 83

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

4) Représenter graphiquement cette erreur en fonction de n+1 et des valeurs du paramètre r.

Voici le script Matlabr

clear␣all␣;␣clc␣;␣close␣all␣;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%␣19/11/2019␣Samir␣KENOUCHE␣:␣ALGORITHME␣PERMETTANT
%␣L’IMPLEMENTATION,␣SOUS␣MATLAB,␣DE␣LA␣METHODE␣DES␣DIFFERENCES␣FINIES␣EN
%␣DIMENSION␣1

fx␣=␣@(x,r)␣(r.ˆ2.*exp(r.*x).*(x-1))./(1␣-␣exp(r))␣+␣r.*x␣;
n␣=␣64␣;␣r␣=␣1/2␣;␣h␣=␣1/(n+1)␣;␣xh␣=␣0␣:␣h:␣1␣;
cx␣=␣inline(’x’)␣;␣funex␣=␣@(x,r)␣x-(1-exp(r.*x))./(1-␣exp(r))␣;

fn␣=␣h.ˆ2.*fx(xh(1:n),r)␣;␣cx␣=␣cx(xh(1:n)).*r␣;

sur_diag␣=␣diag(ones(n-1␣,1)␣,1)␣;␣sur_diag(sur_diag␣==␣1)␣=␣-1␣;
des_diag␣=␣diag(ones(n-1␣,1)␣,␣-1)␣;␣des_diag(des_diag␣==␣1)␣=␣-1␣;
in_diag␣=␣diag(ones(n␣,1))␣;␣in_diag(in_diag␣==␣1)␣=␣2␣;

A1␣=␣sur_diag␣+␣des_diag␣+␣in_diag␣;␣a1␣=␣diag(ones(n␣,1))␣;
a1(a1␣==␣1)␣=␣cx␣;

sur_diag_a2␣=␣diag(ones(n-1␣,1)␣,1)␣;␣sur_diag_a2(sur_diag_a2␣==␣1)␣=␣1␣;
des_diag_a2␣=␣diag(ones(n-1␣,1)␣,␣-1)␣;␣des_diag_a2(des_diag_a2␣==␣1)␣=␣-1␣;
in_diag_a2␣=␣diag(ones(n␣,1))␣;␣in_diag_a2(in_diag_a2␣==␣1)␣=␣0␣;
Cours complet est disponible sur mon site web : [Link]

a2␣=␣sur_diag_a2␣+␣des_diag_a2␣+␣in_diag_a2␣;
A2␣=␣(r*h/2).*(a1*a2)␣;␣An␣=␣A1␣+␣A2␣;

un␣=␣fn*inv(An)␣;␣␣%␣RESOLUTION␣DU␣SYSTEME␣LINEAIRE
uh␣=␣[0␣,␣un,␣0]␣;␣%␣SOL.␣FINALE␣-␣PRISE␣EN␣COMPTE␣DES␣CONDITIONS␣INITIALES

%%%%%%%%%%%%%%%%%%%%%%%%%%%%␣AFFICHAGE␣GRAPHIQUE␣%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(’color’,[1␣1␣1])

plot(xh,uh,’o’,’MarkerSize’,7,’LineWidth’,1)␣;␣hold␣on␣;
xk␣=␣0:0.001:1␣;␣plot(xk,funex(xk,r),’r’,’LineWidth’,1.2)
axis([-0.1␣1.1␣-0.005␣0.07])␣;

ih␣=legend(’SOLUTION␣NUMERIQUE’,’SOLUTION␣EXACTE’)␣;
set(ih,’Interpreter’,’none’,’Location’,’South’,’Box’,’on’,...
␣␣␣␣’Color’,’none’)␣;␣xlabel(’x’,’FontSize’,12)␣;␣ylabel(’u(x)’,’FontSize’,12)

msg1␣=␣strcat(’r=␣’,␣num2str(r))␣;
gtext(msg1)␣%␣cliquer␣sur␣la␣figure␣pour␣afficher␣:␣msg1
msg2␣=␣strcat(’n=␣’,␣num2str(n))␣;
gtext(msg2)␣%␣cliquer␣sur␣la␣figure␣pour␣afficher␣:␣msg2

Chapitre3 Équations aux dérivées partielles et fonctions de Green 84

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Ci-dessous, la visualisation graphique de la solution obtenue.

r = 0.5
n = 64

u(x)

Figure 4: Solutions exacte et numérique obtenues par différences finies

clear␣all␣;␣clc␣;␣close␣all␣;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%␣␣18/11/2019␣Samir␣KENOUCHE␣:␣ALGORITHME␣PERMETTANT
%␣L’ANALYSE␣DE␣L’ERREUR␣EN␣FONCTION␣DU␣NOMBRES␣DE␣SOUS-INTERVALLES␣DE
%␣DISCRETISATION
fx␣=␣@(x,r)␣(r.ˆ2.*exp(r.*x).*(x-1))./(1␣-␣exp(r))␣+␣r.*x␣;␣␣r␣=␣1/2␣;
Cours complet est disponible sur mon site web : [Link]

cx1␣=␣inline(’x’)␣;␣funex␣=␣@(x,r)␣x-(1-exp(r.*x))./(1-␣exp(r))␣;
n␣=␣5␣:␣5␣:␣120␣;

for␣ik␣=␣1:length(n)

␣␣h␣=␣1/(n(ik)+1)␣;␣xh␣=␣0:h:␣1␣;
␣␣fn␣=␣h.ˆ2.*fx(xh(1:n(ik)),r)␣;␣cx␣=␣cx1(xh(1:n(ik))).*r␣;

␣sur_diag␣=␣diag(ones(n(ik)-1␣,1)␣,1)␣;␣sur_diag(sur_diag␣==␣1)␣=␣-1␣;
␣des_diag␣=␣diag(ones(n(ik)-1␣,1)␣,␣-1)␣;␣des_diag(des_diag␣==␣1)␣=␣-1␣;
␣in_diag␣=␣diag(ones(n(ik),␣1))␣;␣in_diag(in_diag␣==␣1)␣=␣2␣;

␣A1␣=␣sur_diag␣+␣des_diag␣+␣in_diag␣;␣a1␣=␣diag(ones(n(ik)␣,1))␣;
␣a1(a1␣==␣1)␣=␣cx␣;

␣sur_diag_a2␣=␣diag(ones(n(ik)-1␣,1)␣,1)␣;␣sur_diag_a2(sur_diag_a2␣==␣1)␣=␣1␣;
␣des_diag_a2␣=␣diag(ones(n(ik)-1␣,1)␣,␣-1)␣;
␣des_diag_a2(des_diag_a2␣==␣1)␣=␣-1␣;
␣in_diag_a2␣=␣diag(ones(n(ik)␣,1))␣;␣in_diag_a2(in_diag_a2␣==␣1)␣=␣0␣;

␣a2␣=␣sur_diag_a2␣+␣des_diag_a2␣+␣in_diag_a2␣;
␣A2␣=␣(r*h/2).*(a1*a2)␣;␣An␣=␣A1␣+␣A2␣;

␣un␣=␣fn*inv(An)␣;␣%␣RESOLUTION␣DU␣SYSTEME␣LINEAIRE

Chapitre3 Équations aux dérivées partielles et fonctions de Green 85

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

␣uh␣=␣[0␣,␣un,␣0]␣;␣%␣SOL.␣FINALE␣-␣PRISE␣EN␣COMPTE␣DES␣CONDITIONS␣INITIALES

␣err(ik)␣=␣max(abs(uh␣-␣funex(xh,r)))␣;

end

figure(’color’,[1␣1␣1])␣;␣hold␣on␣;␣box␣on␣;
loglog(n+1,err,’o’,’MarkerSize’,7,’LineWidth’,1)␣;␣%␣ECHELLE␣LOGARITHMIQUE
xlabel(’n␣+␣1’,’FontSize’,12)␣;␣ylabel(’Erreur␣absolue’,’FontSize’,12)␣;

Figure 5: Erreur absolue versus le nombre d’intervalles de discrétisation.

Cours complet est disponible sur mon site web : [Link]

clear␣all␣;␣clc␣;␣close␣all␣;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%␣␣18/11/2019␣Samir␣KENOUCHE␣:␣ALGORITHME␣PERMETTANT
%␣L’ANALYSE␣DE␣L’ERREUR␣EN␣FONCTION␣DU␣NOMBRES␣DE␣SOUS-INTERVALLES␣DE
%␣DISCRETISATION
fx␣=␣@(x,r)␣(r.ˆ2.*exp(r.*x).*(x-1))./(1␣-␣exp(r))␣+␣r.*x␣;
cx1␣=␣inline(’x’)␣;␣funex␣=␣@(x,r)␣x-(1-exp(r.*x))./(1-␣exp(r))␣;
n␣=␣64␣;

for␣r␣=␣1:40

␣␣h␣=␣1/(n␣+␣1)␣;␣xh␣=␣0:h:␣1␣;
␣␣fn␣=␣(h.ˆ2).*fx(xh(1:n),r)␣;␣cx␣=␣cx1(xh(1:n)).*r␣;

␣sur_diag␣=␣diag(ones(n␣-␣1␣,1)␣,1)␣;␣sur_diag(sur_diag␣==␣1)␣=␣-1␣;
␣des_diag␣=␣diag(ones(n␣-␣1␣,1)␣,␣-1)␣;␣des_diag(des_diag␣==␣1)␣=␣-1␣;
␣in_diag␣=␣diag(ones(n,␣1))␣;␣in_diag(in_diag␣==␣1)␣=␣2␣;

␣A1␣=␣sur_diag␣+␣des_diag␣+␣in_diag␣;␣a1␣=␣diag(ones(n,␣1))␣;
␣a1(a1␣==␣1)␣=␣cx␣;

␣sur_diag_a2␣=␣diag(ones(n␣-␣1␣,1)␣,1)␣;␣sur_diag_a2(sur_diag_a2␣==␣1)␣=␣1␣;
␣des_diag_a2␣=␣diag(ones(n␣-␣1␣,1)␣,␣-1)␣;

Chapitre3 Équations aux dérivées partielles et fonctions de Green 86

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

␣des_diag_a2(des_diag_a2␣==␣1)␣=␣-1␣;
␣in_diag_a2␣=␣diag(ones(n␣,1))␣;␣in_diag_a2(in_diag_a2␣==␣1)␣=␣0␣;

␣a2␣=␣sur_diag_a2␣+␣des_diag_a2␣+␣in_diag_a2␣;
␣A2␣=␣(r*h/2).*(a1*a2)␣;␣An␣=␣A1␣+␣A2␣;

␣un␣=␣fn*inv(An)␣;␣%␣RESOLUTION␣DU␣SYSTEME␣LINEAIRE
␣uh␣=␣[0␣,␣un,␣0]␣;␣%␣SOL.␣FINALE␣-␣PRISE␣EN␣COMPTE␣DES␣CONDITIONS␣INITIALES

␣err_r␣=␣max(abs(uh␣-␣funex(xh,␣r)))␣;

␣figure(1)␣;␣hold␣on␣;␣box␣on␣;
␣semilogy(r,err_r,’o’,’MarkerSize’,7,’LineWidth’,1)
␣%␣ECHELLE␣SEMI-LOGARITHMIQUE
␣xlabel(’Valeur␣de␣r’,’FontSize’,12)␣;
␣ylabel(’Erreur␣absolue’,’FontSize’,12)␣;

end
Cours complet est disponible sur mon site web : [Link]

Figure 6: Erreur absolue en fonction du paramètre r.

Supplément : avec une démarche analogue on peut résoudre également le problème de la flexion simple
dont la formulation mathématique est donnée par :

( 00
−u (x) + c(x) u(x) = f (x), x ∈ [a, b]
P: (176)
u(a) = α et u(b) = β

En utilisant la formule des différences finies centrées, le problème devient :


u(xi−1 ) − 2 u(xi ) + u(xi+1 )
− + c(xi ) ui = f (xi ) i ∈ {1, ..., n}

h2 (177)

u0 = α et un = β

Chapitre3 Équations aux dérivées partielles et fonctions de Green 87

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

En adoptant une notation matricielle :

Ah uh = bh (178)
Avec,
 
2 −1 0 ... 0
c(x1 ) 0... 0
 
. .. 
−1 . .

−1 2 . 
. .. 
c(x2 ) . .
  
(0)  0 .  (0) 1  .. .. .. 
Ah = Ah + 
 . .. ..
 et Ah = 2 0 . . . 0
 (179)
 .. . .

0  h  .
..

 .
 . . −1 2 −1

0 ... 0 c(xn )
0 ... 0 −1 2

Les deux matrices ci-dessus peuvent se combiner pour donner :

 
2 + c(x1 ) h2 −1 0 ... 0
 .. .. 

 −1 2 + c(x2 ) h2 −1 . 
 .
1  .. .. .. 
Ah = 2  0 . . 0 .  (180)
h  
 .. ..

. . −1 2 + c(xn−1 ) h2 −1
 
 
0 ... 0 −1 2 + c(xn ) h2

f (x1 ) + α h−2
 

 f (x2 ) 


bh = 
 .
..

 (181)
Cours complet est disponible sur mon site web : [Link]


 f (xn−1 ) 
 

f (xn ) + β h−2

Le vecteur des solution (inconnues) aux points xi est donné par :

 
u1



 u2

uh = 
 ..
 ∈ Rn (182)
.


un−1 
 

un

2) Problème non-linéaire: considérons le problème non-linéaire suivant :

00




−u (x) + x u(x)3 = f (x), 0<x<L

P: (183)



u(0) = α
 u(L) = β

Contrairement au cas linéaire (problème (153)) où le terme x u(x)3 n’existe pas, pour le problème (183) nous
avons une relation non-linéaire entre la source de la déformation f (x) et l’amplitude de la déformation u(x).
Autrement dit, si j’applique par exemple une force f (x) deux fois plus grande, l’amplitude de la déformation

Chapitre3 Équations aux dérivées partielles et fonctions de Green 88

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

u(x) n’est pas doublée. En appliquant la formule des différences finies centrées il vient :

ui−1 − 2 ui + ui+1

 −


2
+ xi u3i = f (xi ) i ∈ {1, ..., n}
h (184)



u0 = α et un+1 = β
Contrairement au problème (153), ici nous cherchons à résoudre un système de n équations non-linaires.
Il existe dans la littérature plusieurs méthodes itératives pour effectuer ce calcul. On citera les méthodes
de Newton, de fausse position ou de la sécante. Je renvoie, les lecteurs intéressés par ces méthodes, à mon
cours d’analyse numérique que je dispense aux deuxièmes années des filières physique et chimie. Ce cours est
disponible en version pdf. Pour n = 4, on obtient le système d’équations :

−u0 + 2 u1 − u2



2
+ x1 u31 = f (x1 ), i=1
h

−α + 2 u1 − u2 + h2 x1 u31 − h2 f (x1 ) = 0

 

 

 


 −u1 + 2 u2 − u3 

+ x2 u32 = f (x2 ), i=2
 
2 3 2
 −u1 + 2 u2 − u3 + h x2 u2 − h f (x2 ) = 0

 


 h2 
⇒
 −u2 + 2 u3 − u4 
−u2 + 2 u3 − u4 + h2 x3 u33 − h2 f (x3 ) = 0
+ x3 u33 = f (x3 ),
 
i=3

 

2
 



 h 



−u3 + 2 u4 − β + h2 x4 u34 − h2 f (x4 ) = 0

 
−u3 + 2 u4 − u5


+ x4 u34 = f (x4 ),



2
i=4
h

En posant u0 = 0 et un+1 = 0 (pas de déformations aux extrémités) on obtient :

2 u1 − u2 + h2 x1 u31 − h2 f (x1 ) = 0

Cours complet est disponible sur mon site web : [Link]







2 3 2
 −u1 + 2 u2 − u3 + h x2 u2 − h f (x2 ) = 0




−u2 + 2 u3 − u4 + h2 x3 u33 − h2 f (x3 ) = 0










−u3 + 2 u4 + h2 x4 u34 − h2 f (x4 ) = 0


On peut par exemple résoudre ce système par la méthode de Newton, c’est une méthode itérative d’ordre
deux donc elle converge rapidement. Nous allons illustrer cette méthode à l’aide d’un exemple. Soit à résoudre
le système d’équations non-linéaires suivant :

u
 f1 (u1 , u2 ) = 2 u1 − u2 + e 1 = 0

f (U ) = (185)
 f (u , u ) = −u + 2 u + eu2 = 0

2 1 2 1 2

Nous cherchons u1 et u2 telle que f (U ) = 0, avec le vecteur U = (u1 , u2 )T . Le schéma numérique de la

méthode de Newton pour résoudre le système (185) est :

f (U )
Uk+1 = Uk −
∇f (U )
f (u1 , u2 )
[uk+1 k+1 k k
1 , u2 ] = [u1 , u2 ] −
∇f (u1 , u2 )

Chapitre3 Équations aux dérivées partielles et fonctions de Green 89

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Avec ∇f (u1 , u2 ) est la matrice Jacobienne :
∂f1 ∂f1
 
 ∂u ∂u2 
 1 
∇f (u1 , u2 ) = 


 (186)
 ∂f2 ∂f2 
∂u1 ∂u2
On donne une valeur initiale au vecteur Uk ensuite on calcule les successeurs Uk+1 . On cesse les itérations
f (U )
une fois le test d’arrêt est positif pour une tolérance donnée, par exemple | − ∇f (U ) | < .

3) En dimension 2: il existe une myriade de problèmes en physique et en chimie admettant comme formu-
lation mathématique, une équation aux dérivées partielles. Rappelons que cette dernière exprime une relation
fonctionnelle dont l’inconnue est une fonction de plusieurs variables. Dans l’équation même apparait la fonction
de plusieurs variables recherchée ainsi que ses dérivées partielles. Soit ρ : (x, t) ∈ [0, 1]×R+ 7−→ ρ(x, t) ∈ R une
fonction continue, nous considérons un problème parabolique consistant à déterminer u : (x, t) ∈ [0, 1]×R+ 7−→
u(x, t) ∈ R qui satisfait :

∂ 2 u(x, t)

∂u(x, t)
−


 τ c p k = ρ(x, t)
∂t ∂x2





P: (187)

 u(x, 0) = u0 (x) 0≤x≤1






u(0, t) = u(a, t) = 0 0 ≤ t

C’est l’équation de la diffusion de la chaleur, avec ρ(x, t) est la source de chaleur. Les constantes positives
τ , cp et k, caractéristiques du matériau en question, représentent respectivement la densité volumique, la
chaleur spécifique massique et la conductivité thermique. Afin d’alléger les écritures ces coefficients sont pris
Cours complet est disponible sur mon site web : [Link]

égaux à l’unité. A partir de ce problème, nous cherchons à déterminer la quantité de chaleur fournie au point
x à l’instant t. A partir d’un développement de Taylor on démontre ces approximations des dérivées partielles :

∂u(x, t) u(xi , tj+1 ) − u(xi , tj )

≈ + O(∆t) (188)
∂t ∆t
∂u(x, t) u(xi , tj ) − u(xi , tj−1 )
≈ − O(∆t) (189)
∂t ∆t

∂ 2 u(x, t) u(xi−1 , tj ) − 2 u(xi , tj ) + u(xi+1 , tj )

2
≈ + O(∆x2 ) (190)
∂x ∆x2
Nous utiliserons ces schémas des différences finies afin d’approcher u(x, t) du problème ci-dessus. Nous cher-
chons des approximations u(i, j) de la solution exacte u(xi , tj ) aux nœuds (xi , tj ) = (i×∆x, j×∆t)i,j={0,1,...,n+1} .
Cette discrétisation définie un maillage ou une grille selon le domaine Ω = [i × ∆x , j × ∆t]2 . La condition
initiale u(x, 0) ' u(i × ∆x, 0) = 0 signifie que la quantité de chaleur u(i × ∆x, 0) est connue sur chaque nœud
xi = i×∆x à l’instant initial (t = 0). D’un autre côté, les conditions aux limites u(0, j ×∆t) = u(a, j ×∆t) = 0
signifient que la quantité de chaleur, apportée aux limites (ou aux bords) du domaine Ω, est nulle. En effet,
à partir des équations (189), (190) et en posant u(xi , tj ) ' u(i, j) on obtient :

u(i, j) − u(i, j − 1) u(i − 1, j) − 2u(i, j) + u(i + 1, j)

− = ρ(i, j)
∆t ∆x2
∆x2 u(i, j) − ∆x2 u(i, j − 1) − ∆t u(i − 1, j) + 2 ∆t u(i, j) − ∆t u(i + 1, j) = ρ(i, j) ∆t ∆x2

Chapitre3 Équations aux dérivées partielles et fonctions de Green 90

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Tenant compte des conditions aux limites u(0, j × ∆t) = u(a, j × ∆t) = 0 il vient :

∆x2 u(i, j) − ∆x2 u(i, j − 1) − ( −(

1, j) + 2 ∆t u(i, j) − ( 1, j) = ρ(i, j) ∆t ∆x2
( (
∆t(u(i ∆t(u(i +(
( (( (
((

(∆x2 + 2 ∆t) u(i, j) = ∆x2 u(i, j − 1) + ρ(i, j) ∆t ∆x2

2 ∆t
(I + ) u(i, j) = u(i, j − 1) + ρ(i, j) ∆t
∆x2
An ∆t
(I + ) u(i, j) = u(i, j − 1) + ρ(i, j) ∆t
∆x2

Tenant compte de la condition initiale u(x, 0) ' u(i × ∆x, 0) = 0 le schéma numérique final devient :
An ∆t
(I +
) u(i, j + 1) = u(i, j) + ρ(i, j) ∆t
∆x2
Avec I est la matrice identité. La matrice An est tridiagonale, symétrique et définie positive valant :

 
2 −1 0 ... 0
. .. 
−1 . .

−1 2 . 
 
 .. .. ..  n n
An = 
 0 . . . 0 ∈R ×R

 .
..

 .
 . . −1 2 −1

0 ... 0 −1 2

Voici le script Matlabr

Cours complet est disponible sur mon site web : [Link]

clear␣all␣;␣clc␣;␣close␣all␣;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%␣LE␣04/12/2019␣Samir␣KENOUCHE␣:␣RESOLUTION␣DE␣L’EQUATION␣DE␣LA␣CHALEUR␣PAR
%␣DIFFERENCES␣FINIES␣EN␣DIMENSION␣2

nb␣=␣10␣;␣pas_temps␣=␣0.01␣;␣pas_x␣=␣1/(nb+1)␣;␣T␣=␣0.7␣;
xi␣=␣0:␣pas_x␣:1␣;␣ti␣=␣0:␣pas_temps␣:T␣;␣ux␣=␣sin(pi.*xi)␣;␣%␣Condition
␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣%␣Initiale
u␣=␣zeros([numel(xi)-2␣numel(ti)])␣;␣u(:,␣1)␣=␣ux(2:end-1)’␣;
fx␣=␣-2␣+␣6.*xi␣;␣fn␣=␣zeros(size(u))␣;
mat_diag␣=␣2*diag(ones(nb,1))-diag(ones(nb-1␣,1),1)-diag(ones(nb-1,1),-1)␣;
my_mat␣=␣(eye(nb)+(pas_temps/pas_x.ˆ2).*mat_diag)␣;␣it␣=␣1␣;

while␣␣it␣<␣numel(ti)

␣␣␣␣␣fn(:,␣it)␣=␣fx(2:end-1)’␣;
␣␣␣␣␣u(:,␣it+1)␣=␣my_mat␣\(u(:,it)␣+␣pas_temps*fn(:,it))␣;␣%␣DIVISION␣A␣GAUCHE
␣␣␣␣␣it␣=␣it␣+␣1␣;
end

ui␣=␣cat(1,zeros([1␣numel(ti)]),u,zeros([1␣numel(ti)]));␣%␣AJOUT␣DES␣CONDITION
␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣%␣INITIALE␣ET␣AU␣LIMITE
fig1␣=␣figure(’color’,[1␣1␣1])␣;␣[xn,␣tn]␣=␣meshgrid(xi,␣ti)␣;

Chapitre3 Équations aux dérivées partielles et fonctions de Green 91

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

surf(xn,␣tn,␣ui’)␣;␣xlabel(’Distance’)␣;␣ylabel(’Temps’)␣;␣view(114,␣18)␣;

fig2␣=␣figure(’color’,[1␣1␣1])␣;␣contourf(xn,␣tn,␣ui’)␣;
xlabel(’Distance’)␣;␣ylabel(’Temps’)␣;
%%%%%%%%%%%%%%%%%%%%%%%%%%%␣PARTIE␣:␣ANIMATION␣%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(’Renderer’,’zbuffer’)␣;␣[xn,␣tn]␣=␣meshgrid(xi,␣ti)␣;
colormap(jet)␣;␣xlabel(’Distance’)␣;␣ylabel(’Temps’)␣;

for␣in␣=␣1:30
␣␣␣␣surf(exp(-0.008*in)*ui’.ˆ2,ui’)
␣␣␣␣my_animation(in)␣=␣getframe␣;
end

0.8

0.6 Temps

0.4

0.2
0

0 0.5 e
0 0.1 0.2 a nc
Cours complet est disponible sur mon site web : [Link]

0.3 0.4 0.5 0.6 st

Temps 0.7 1
Di

Figure 7: Surface de la solution approchée

0.6

0.5

0.4
Temps

0.3

0.2

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Distance
Figure 8: Courbe de niveaux de la solution approchée

Chapitre3 Équations aux dérivées partielles et fonctions de Green 92

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

III. Fonctions de Green
Considérons un système physique ou chimique, par exemple un oscillateur mécanique (peut être une liaison
chimique) oscillant autour d’une position d’équilibre. Les oscillations d’amplitude u(t) (sortie, ou réponse du
système) sont causées par une excitation φ(t) (ou source des vibrations de la liaison chimique). D’un point
de vue mathématique, la réponse u(t) est une fonctionnelle de la source φ(t) qui s’écrit donc sous la forme :
u[φ(t)]. Nous rappelons que la différentielle d’une fonctionnelle quelconque F[f (t)] s’écrit comme suit :
Z b δF[f (t)]
δF[f (t)] = F[f (t) + δf (t)] − F[f (t)] = δf (t) dt (191)
a δf (t)
Appliquons désormais cette formule a notre fonctionnelle u[φ(t)], il vient :
Z b δu[φ(t)]
δu[φ(t)] = δφ(τ ) dτ, ∀t > τ (192)
a δφ(τ )
| {z }
h0 (t,τ )

L’introduction de la variable temporelle τ est justifiée par le principe de causalité 9 , qui stipule que la cause
φ(τ ) est systématiquement antérieur à l’effet ou à la réponse u(t), de sorte que nous avons toujours t > τ .
Cela traduit une relation chronologique entre l’excitation et la réponse du système en question.
δu[φ(t)]
En revanche, la quantité h0 (t, τ ) = exprime la réponse propre ou intrinsèque du système considéré.
δφ(τ )
La connaissance de cette quantité mathématique est fondamentale, car elle permet de quantifier la réponse (ou
la sortie) u(t) quelque soit l’excitation (ou la source) φ(t). Cela nous amène à écrire le produit de convolution 10
suivant :
Z b
u(t) = φ(τ ) h0 (t, τ ) dτ (193)
a

Dans ce cas de figure, l’excitation φ(t) est écrite sous forme d’une superposition (ou une somme) d’impulsions
Cours complet est disponible sur mon site web : [Link]

de Dirac, selon :
Z b
φ(t) = φ(τ ) δ(t, τ ) dτ (194)
a

Il faut noter pour que l’équation (192) soit justifiée, le système considéré doit être stable. En effet, un système
donné est dit stable si, en lui appliquant une excitation bornée quelconque, la réponse reste bornée par une
constante, noté m. Traduit en terme mathématique, cela implique :
Z b
h0 (t, τ ) dτ < m (195)
a
Une autre notion fondamentale est l’invariance par translation dans le temps. La plupart des systèmes
physiques respecte cette invariance, stipulant que si l’on retarde l’excitation de δτ , alors la réponse est
également retardée de δτ . Cela se traduit par l’écriture :

h0 (t, τ ) = h0 (t − τ ) (196)

9. La causalité constitue une contrainte majeure et inviolable pour la formalisation de nombreux problèmes en physique et
en chimie. Il es résulte ce qui suit : t < τ ⇒ h0 (t, τ ) = 0
10. Cette relation est linéaire, car elle vérifie les propriétés suivantes :
1) ∀α ∈ R : α φ(t) ⇒ α u(t).
2) ∀α, β ∈ R : α φ1 (t) + β φ2 (t) ⇒ α u1 (t) + β u2 (t).

Chapitre3 Équations aux dérivées partielles et fonctions de Green 93

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

En effet, pour une excitation décalée φ2 (t − δτ ) = φ1 (t) et d’après ce que nous avons vu précédemment,
nous pouvons écrire :
Z b Z b
φ1 (τ ) h0 (t, τ ) dτ = φ2 (τ − δτ ) h0 (t, τ ) dτ (197)
a a
0
En considérant (196) et en opérant le changement de variable suivant t = τ − δτ , il vient :
Z b Z b 0 0
φ1 (τ ) h0 (t, τ ) dτ = φ2 (t ) h0 (t − τ ) dt ⇔ u2 (t) = u1 (t − δτ ) (198)
a a
Cela signifie que le système ne change pas ses caractéristiques, à mesure que le temps passe. Compte tenue de
toutes les propriétés (causalité, linéarité et invariance) énumérées précédemment, la réponse se réécrit selon :
Z b
u(t) = φ(τ ) h0 (t − τ ) dτ (199)
a

L’évolution temporelle associée à la grandeur u(t) est ainsi simplement proportionnelle à la fonction h0 (t−τ ).
Une fois ces propriétés fondamentales sont rappelées, nous présenterons dans ce qui suit, le cadre théorique
des fonctions de Green. Un problème aux limites d’ordre n est formulé mathématiquement comme suit :

 L̂(t) u(t) = φ(t),
 t ∈ [a, b]
(200)

 U (u) = γ ,
j j j = 1, m
Où,
n−1
αij u(i) (a) + βij u(i) (b), ∀αij , βij ∈ R,
X
Uj (u) = n≤m (201)
i=0

Avec L̂ est un opérateur de dérivation linéaire. L’équation (201) traduit les conditions aux limites imposées
Cours complet est disponible sur mon site web : [Link]

à la solution u(t). Le problème homogène associé à (200) s’écrit selon :


 L̂(t) u(t) = 0,
 t ∈ [a, b]
(202)

 U (u) = 0,
j j = 1, m
Théorème de Fredholm : Le problème aux limites (200) admet une solution unique si et seulement si, le
problème homogène L̂(t) u(t) = 0, t ∈ [a, b], Uj (u) = 0, admet comme solution triviale u(t) = 0.

Définition : G(t, τ ) est une fonction de Green du problème homogène, alors ∃ une fonction continue φ(t)
telle que :
Z b
u(t) = G(t, τ ) φ(τ ) dτ
a

soit une solution particulière du problème aux limites (200). De plus, la fonction G(t, τ ) doit satisfaire les
propriétés suivantes :

1) Pour chaque τ ∈ [a, b], la fonction t 7−→ G(t, τ ) est solution de l’équation homogène L̂(t) u(t) =
0, Uj (u) = 0, ∀t ∈ [a, τ ] et t ∈ [τ, b]. Cela se traduit mathématiquement par :

∂ n G(t, τ ) ∂ n−1 G(t, τ ) ∂G(t, τ )

an (t) n
+ a n−1 (t) n−1
+ · · · + a1 (t) + a0 (t) G(t, τ ) = 0 (203)
∂t ∂t ∂t

Chapitre3 Équations aux dérivées partielles et fonctions de Green 94

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

2) Pour chaque τ ∈ [a, b], la fonction t 7−→ G(t, τ ) vérifie les conditions aux limites homogènes Uj (G(t, τ )) =
0, j = 1, m. Autrement dit,
n−1
" #
∂ i G(a, τ ) i
j ∂ G(b, τ )
αij ∀αij , βij ∈ R,
X
+ βi = 0, n≤m (204)
i=0
∂ti ∂ti
Cela se démontre en considérant (201), selon :

0 0
U1 = α01 u(a) + β01 u(b) + α11 u (a) + β11 u (b)
Z b h 0 0
i
= α01 G(a, τ ) + β01 G(b, τ ) + α11 G (a, τ ) + β11 G (b, τ ) φ(τ ) dτ
a
Z b
= U1 (G(t, τ )) φ(τ ) dτ = 0
a | {z }
car U1 (G(t,τ ))=0

3) La fonction de Green est continue à t = τ :

lim G(t, τ ) = lim G(t, τ ) ⇔ lim G(t, τ ) = lim G(t, τ ) (205)

t→τ + t→τ −  τ >t τ <t
| {z } | {z }
τ+ τ−

4) Une propriété intéressante de la fonction de Green, est la présence d’une discontinuité d’une amplitude
1
valant de sa dérivée d’ordre (n − 1), à t = τ :
an (t)
∂Gn−1 (t, τ ) ∂Gn−1 (t, τ ) 1
lim+ t>τn−1 − lim− t<τn−1 = (206)
t→τ ∂t t→τ ∂t an (t)
Avec an (t) étant le coefficient de la dérivée d’ordre n de u(t).
Cours complet est disponible sur mon site web : [Link]

Les fonctions de Green permettent de transformer une équation différentielle (ou aux dérivées partielles) en
une équation intégrale, généralement plus simple à résoudre. Ces fonctions sont largement utilisées en Chimie
et Physique quantique. L’équation (200) est explicitée sous la forme suivante :
n n
X
(k)
X dk
ak (t) u (t) = φ(t) ⇐⇒ L̂(t) u(t) = φ(t), L̂(t) = ak (t) (207)
k=0 k=0
dtk
Les coefficients ak (t) peuvent aussi êtres constants. La fonction de Green associée à l’équation ci-dessus, est
obtenue en remplaçant la source φ(t) par une impulsion appliquée à un instant τ > t, soit δ(t − τ ). La source
φ(t) pouvant être une force agissant sur une particule par exemple, ou encore une source de chaleur ou de
vibration ... etc. Quelque soit la nature de cette source, elle peut être vue comme une somme d’impulsions
appliquées à différents instants τ > t. Cela se traduit par mathématiquement par l’écriture :
Z ∞
φ(t) = φ(τ ) δ(t − τ )dτ (208)
0

Où φ(τ ) est la fonction de poids associée à chaque impulsion. De cette façon, la source φ(t) est totalement
reconstituée par le second membre de l’équation (208). Revenons maintenant à notre équation différentielle
(207), si l’on note G(t, τ ) la fonction de Green 11 alors :
n
X ∂ k G(t, τ )
ak (t) = δ(t − τ ) ⇐⇒ L̂(t) G(t, τ ) = δ(t − τ ) (209)
k=0
∂tk

11. Avec τ est une variable factice, ou dummy variable en anglais

Chapitre3 Équations aux dérivées partielles et fonctions de Green 95

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

En développant (209) à l’ordre deux et en optant pour des coefficients constants, il vient :
00 0
a2 G (t, τ ) + a1 G (t, τ ) + a0 G(t, τ ) = δ(t − τ ) (210)
Résoudre l’équation (210), revient à déterminer d’abord la solution générale (L̂(t) G(t, τ ) = 0), à laquelle
nous ajoutons la solution particulière (L̂(t) G(t, τ ) = δ(t − τ )). Multiplions (210) par la fonction de poids φ(τ ),
il vient :

∂2 ∂
2
[φ(τ ) G(t, τ )] + a1
a2 [φ(τ ) G(t, τ )] + a0 [φ(τ ) G(t, τ )] = φ(τ ) δ(t − τ ) (211)
∂t ∂t
Intégrons les deux membres de cette équation et nous obtenons :

∂2 ∞ ∂ ∞ ∞ ∞
Z Z Z Z
a2 2 [φ(τ ) G(t, τ )] dτ + a1 [φ(τ ) G(t, τ )] dτ + a0 [φ(τ ) G(t, τ )] dτ = φ(τ ) δ(t − τ )dτ
∂t 0 ∂t 0 0
|0 {z }
φ(t)

De cette manière nous avons totalement satisfait l’équation différentielle (207). D’après ce dernier résultat,
nous comprenons que l’intégrale :
Z ∞
up (t) = φ(τ ) G(t, τ ) dτ (212)
0

constitue une solution particulière de l’équation différentielle (207). La solution globale (générale + parti-
culière) sera donnée par :
Z ∞
u(t) = uh (t) + up (t) = uh (t) + φ(τ ) G(t, τ ) dτ (213)
0

Où uh (t) dénote la solution générale de l’équation homogène (sans second membre, L̂(t) u(t) = 0) associée.
Ainsi, si l’on connait la fonction de Green G(t, τ ) d’une équation différentielle, une solution particulière s’obtient
Cours complet est disponible sur mon site web : [Link]

par l’équation intégrale up (t). Toute la problématique consiste donc à déterminer la ”bonne” fonction de Green
du phénomène physique ou chimique étudié. Toutefois, la détermination de G(t, τ ) est totalement tributaire
de la connaissance de la solution de l’équation homogène uh (t).

Par ailleurs, nous aurions pu démontrer l’équation (213), en utilisant une des propriétés de la fonction de
Dirac à savoir :
Z ∞
φ(τ ) δ(t − τ ) dτ = φ(t) (214)
0

D’après (209), nous avons :

Z ∞ Z ∞
L̂(t) G(t, τ ) = L̂(t) φ(τ ) G(t, τ ) dτ = φ(τ ) δ(t − τ ) dτ (215)
0 0
| {z }
φ(t)

L’opérateur différentiel L̂(t) est indépendant de τ , par conséquent et comme précédemment une solution de
(207) s’obtient selon :
Z ∞
up (t) = φ(τ ) G(t, τ ) dτ (216)
0

Il est important de rappeler que la fonction G(t, τ ) satisfait pleinement les conditions aux limites imposées
à la solution u(t).

Exercice d’application : on se propose d’étudier un problème aux limites du second ordre dont la
formulation mathématique est donnée comme suit :

Chapitre3 Équations aux dérivées partielles et fonctions de Green 96

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

 00
 u (t) = φ(t)
 avec 0 ≤ t ≤ 2π
P:
 u(0) = u(2π) = 1


D’après ce que nous avons discuté précédemment, la solution de (P) s’écrira comme :
Z ∞
u(t) = uh (t) + φ(τ ) G(t, τ ) dτ (217)
0

De plus, à partir de (P), nous écrivons l’équation satisfaite par la fonction de Green associée soit :
00
G (t, τ ) = δ(t − τ ) (218)
Commençons d’abord par résoudre l’équation homogène associée :




a1 t + b1 0≤t<τ
00

G (t, τ ) = 0 ⇒ G(t, τ ) = (219)

 a2 t + b2 τ < t ≤ 2π



Nous allons étudier le cas où t < τ , ou de façon totalement équivalente quand 0 ≤ t < τ d’après la condition
initiale. En effet, tenant compte de la condition u(0) = 1, exprimée dans (P), il vient :

Gt<τ (t, τ ) = a1 t + b1 ⇒ Gt<τ (t = 0, τ ) = b1 = 1 (220)

⇒ Gt<τ (t, τ ) = a1 t + 1 (221)

Étudions également le cas où τ < t, ou de façon totalement équivalente quand τ < t ≤ 2π d’après les
conditions aux limites. Ainsi, tenant compte de la condition u(2π) = 1, exprimée dans (P), il vient :
Cours complet est disponible sur mon site web : [Link]

Gt>τ (t, τ ) = a2 t + b2 ⇒ Gt>τ (t = 2π, τ ) = 2 a2 π + b2 = 1 ⇒ b2 = 1 − 2 a2 π (222)

⇒ Gt>τ (t, τ ) = a2 (t − 2π) − 1 (223)

Nous obtenons la fonction de Green suivante :




a1 t + 1 0≤t<τ

G(t, τ ) = (224)

 a2 (t − 2π) + 1 τ < t ≤ 2π



Étudions désormais le cas où t = τ , nous savons que dans ce cas de figure la fonction G(t, τ ) est continue,
soit :

lim G(t, τ ) = lim G(t, τ ) (225)

τ >t τ <t

Alors,

t=τ ⇒ Gt<τ (t, τ ) = Gt>τ (t, τ ) ⇔ a1 τ + 1 = a2 (τ − 2π) + 1 (226)

Pour déterminer les constantes a1 et a2 , exploitons la discontinuité de G(t, τ ) par rapport à sa dérivée
d’ordre un, soit :
τ + ∂ 2 G(t, τ ) τ + ∂Gt>τ (t, τ ) ∂Gt<τ (t, τ )
Z Z
dt = δ(t − τ ) dt = lim+ − lim− = a2 − a1 = 1 (227)
τ − ∂t2 τ − t→τ ∂t t→τ ∂t

Chapitre3 Équations aux dérivées partielles et fonctions de Green 97

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

⇒ a1 = a2 − 1 (228)
D’après (226) et (228), il en ressort :
τ τ − 2π
(a2 − 1) τ = a2 (τ − 2π) ⇒
et donc a1 = a2 = (229)
2π 2π
Finalement, la fonction de Green associée à P, s’écrit sous la forme :
τ − 2π


 t pour t < τ
2π


G(t, τ ) = (230)

 τ

 (t − 2π)
pour t > τ
2π
Les deux formules apparaissant dans (230) sont les expressions d’une seule et même fonction à savoir
G(t, τ ). Chaque expression est valide dans une région donnée, t < τ ou t > τ . Toutefois, la fonction de
τ
Green G(t, τ ) = (t − 2π) est celle qui répond le plus aux considérations de la physique, sachant que la
2π
réponse du système u(t) intervient après une impulsion appliquée à un instant antérieur. Quelque soit la forme
mathématique de la fonction source φ(t) (qui est connue), la solution finale de notre équation différentielle
sera obtenue par (217) :
(t − 2π)
Z 2π
u(t) = uh + φ(τ ) τ dτ (231)
2π 0 | {z }
fonction donnée
00
Cherchons la solution uh de l’équation homogène u = 0. L’équation caractéristique s’écrit :

λ20 = 0 ⇒ λ0 = 0 (232)
Le discriminant de l’équation caractéristique ∆ = b2 − 4 a c = 0, c’est une racine double et la solution de
Cours complet est disponible sur mon site web : [Link]

l’équation homogène est donnée par :

uh (t) = (c1 + c2 t) eλ0 t ∀c1 , c2 ∈ R (233)

Il es ressort la solution globale suivante :
(t − 2π)
Z 2π
u(t) = (c1 + c2 t) + φ(τ ) τ dτ, ∀c1 , c2 ∈ R (234)
2π 0 | {z }
fonction donnée

Par exemple, si la source φ(t) = sin(w0 t), alors :

(t − 2π)
Z 2π
u(t) = (c1 + c2 t) + sin(w0 τ ) τ dτ, ∀c1 , c2 ∈ R (235)
2π 0
Un intégration par partie nous donne :
Z Z
u dv = u v − v du (236)

Posons,

u=τ ⇒ du = dτ

− cos(w0 t)
dv = sin(w0 τ ) dτ ⇒ v=
w0
Z 2π sin(2π w0 ) − 2π w0 cos(w0 2π)
sin(w0 τ ) τ dτ = (237)
0 w02

Chapitre3 Équations aux dérivées partielles et fonctions de Green 98

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

La solution finale est :
(t − 2π) sin(2π w0 ) − 2π w0 cos(w0 2π)
 
u(t) = (c1 + c2 t) + , ∀c1 , c2 ∈ R (238)
2π w02

A. Fonction de Green de l’équation de Schrödinger

En mécanique quantique, les particules élémentaires, apparaissent pour un observateur donné, se comporter
comme des ondes appelées ondes de Broglie. La mécanique de ces objets devient celle du mouvement des ondes
et les fonctions d’onde sont utilisées pour décrire le comportement des ces systèmes quantiques. Le module
au carré de la fonction d’onde traduit la probabilité qu’une particule existe en un point donné de l’espace.
C’est pourquoi, les fonctions d’ondes sont parfois appelées ondes de probabilité pouvant êtres dispersées par un
potentiel atomique ou nucléaire V (r). Nous cherchons une solution de l’équation de Schrödinger stationnaire,
par exemple, pour le mouvement d’électrons d’énergie E dans un champ d’énergie potentielle V (r) 12 . Si le
potentiel est un diffuseur élastique et les ondes de de Broglie décrivent des particules non relativistes, l’équation
aux dérivées partielles (indépendante du temps) qui décrit le mieux cet effet (de diffusion) est donnée par :

~2 2
∇ ψ(r) + V (r) ψ(r) = E ψ(r)
− (239)
2m
Où k est le nombre d’onde et V (r) est une inhomogénéité qui est responsable de la diffusion de l’onde
ψ(r) qui est donc parfois appelé diffuseur. Cette équation est connue sous le nom d’équation de Schrödinger
d’après le physicien théoricien autrichien Erwin Schrodinger qui l’a postulée dans les années 1920. Après
réarrangement, nous obtenons :
2mE 2m
∇2 ψ(r) + 2
ψ(r) = 2 V (r) ψ(r) (240)
| ~
{z } |~ {z }
k2 Q(r)
h i
⇒ ∇2 + k 2 ψ(r) = Q(r) (241)
Cours complet est disponible sur mon site web : [Link]

La relation (241) a la forme de l’équation de Helmhotz. Notons, cependant, que le terme inhomogène Q(r)
dépend explicitement de la solution ψ(r). Nous commencerons par étudier la solution de la fonction de Green
vérifiant l’équation de Helmholtz inhomogène. Comme pour l’exemple d’application précédent, la fonction de
Green associée à (241) doit être solution de l’équation suivante :
h i
⇒ ∇2 + k 2 G(r) = δ(r) (242)
Ainsi, la solution particulière de (241) est :
Z +∞
ψp (r) = Q(r0 ) G(r − r0 ) dr0 (243)
−∞

Nous pouvons vérifier facilement que ψp (r) soit solution de (242), autrement dit :
h i Z +∞ h i Z +∞
∇2 + k 2 ψp (r) = ∇2 + k 2 G(r − r0 ) Q(r0 )dr0 = δ(r − r0 ) Q(r0 ) dr0 = Q(r) (244)
−∞ −∞

Toute la problématique consiste donc à déterminer la fonction G(r). D’après ce que nous avons déjà vu,
la fonction G(r) pour une équation différentielle donnée représente la réponse à une impulsion δ(r). Nous
allons d’abord résoudre l’équation (242). Pour se faire, nous recourrons à la Transformée de Fourier afin de
transformer l’équation différentielle (242) en une équation algébrique plus simple à résoudre. Multiplions (242)
par le facteur ejs r , il vient :

12. En théorie quantique, il est d’usage d’appeler potentiel tout court, l’énergie potentielle V (r) qui est effectivement le produit
d’un potentiel par une charge.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 99

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Z +∞ h i Z +∞ Z +∞ Z +∞
∇2 + k 2 G(r) ejs r dr = δ(r) ejs r = 1 ⇒ ∇2 G(r) ejs r dr + k 2 G(r) ejs r dr = 1
−∞ −∞ −∞ −∞

Z +∞ Z +∞ Z +∞ Z +∞
⇒ (js)2 G(r) ejs r dr+k 2 G(r) ejs r dr = 1 ⇒ −s2 G(r) ejs r dr+k 2 G(r) ejs r dr = 1
−∞ −∞ −∞ −∞

h iZ +∞
⇒ k 2 − s2 G(r) ejs r dr = 1
−∞
| {z }
g(s)

Finalement, nous obtenons la Transformée de Fourier de la fonction de Green dans l’espace réciproque
1
s ≡ , soit :
r
1
⇒ g(s) = 2 (245)
k − s2
Pour revenir à l’espace réel, nous devons effectuer une Transformée de Fourier inverse. Par définition nous
avons :

1 +∞ 1 +∞ e−js r
Z Z
G(r) = g(s) e−js r ds ⇒ G(r) = ds (246)
(2 π)3 −∞ (2 π)3 −∞ k 2 − s2
La particule évolue dans un espace en trois dimension, pour représenter cette réalité physique nous devons
adopter le système de coordonnées sphériques pour résoudre l’intégrale ci-dessus. L’élément de volume en
coordonnées sphériques s’écrit :

ds = s2 ds sin(θ) dθdφ (247)

Et d’un autre côté, nous avons aussi
Cours complet est disponible sur mon site web : [Link]

e−js r = e−j |s| |r| cos(θ) C’est un produit scalaire entre deux vecteurs (248)
Afin de ne pas alourdir les écritures mathématiques, nous gardons e−js r cos(θ) . Par conséquent l’intégrale
précédente devient :

1 2π +∞ s2 π
Z Z Z
⇒ G(r) = dφ ds sin(θ) e−js r cos(θ) dθ (249)
(2 π)3 0 0 k 2 − s2 0

1 +∞ s2 π
Z Z
⇒ G(r) = − ds sin(θ) e−js r cos(θ) dθ (250)
(2 π)2 0 k 2 − s2 0

Afin de résoudre l’intégrale relative à la variable θ, nous devons procéder au changement de variable suivant :

z = cos(θ) ⇒ dz = − sin(θ) dθ
Évidemment, les bornes d’intégration changent aussi avec la nouvelle variable d’intégration z, selon :

θ=0 ⇒ z = +1
θ=π ⇒ z = −1

1 +∞ s2 1
Z Z
⇒ G(r) = + ds e−js r z dz (251)
(2 π)2 0 k 2 − s2 −1
" #z=1
1 +∞ s2 e−js r z
Z
⇒ G(r) = + ds (252)
(2 π)2 0 k 2 − s2 j sr z=−1

Chapitre3 Équations aux dérivées partielles et fonctions de Green 100

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

h i
1
Z +∞ s ejs r − e−js r
⇒ G(r) = ds (253)
j r (2 π)2 0 k 2 − s2
Par ailleurs, rappelons que :
Z +∞ Z +∞
f (x) dx = 2 f (x) dx (254)
−∞ 0

L’intégrale (253) devient :

h i
j
Z +∞ s ejs r − e−js r
⇒ G(r) = ds (255)
8 r π2 −∞ s2 − k 2
Avec,
1 1
= (256)
(s2 2
−k ) (s − k) (s + k)
"Z #
j +∞ s ejs r +∞ s e−js r
Z
⇒ G(r) = ds − ds (257)
8 r π2 −∞ (s − k) (s + k) −∞ (s − k) (s + k)
Les deux intégrales du second membre sont le contour du demi-cercle à k = ±s. Nous aurons une rotation
dans le sens des aiguilles d’une montre à −k et le chemin inverse à +k. Ces deux intégrales sont résolues en
utilisant la formule de Cauchy suivante :
f (u)
I
= 2π j f (u0 ) (258)
u − u0
En appliquant cette formule, nous obtenons :
" " # " # #
j s ejs r s e−js r
⇒ G(r) = 2π j − (−2π j) (259)
Cours complet est disponible sur mon site web : [Link]

8 r π2 s+k s=k
s−k s=−k

Nous obtenons finalement la fonction de Green pour le problème (242) :

1 jk r
⇒ G(r) = − e (260)
4π r
Ainsi la solution particulière (243), prend la forme 13 :

1 +∞ ejk |r−r0 |
Z
ψp (r) = − Q(r0 ) dr0 (261)
4π r −∞ |r − r0 |
Rappelons que nous avons déjà posé :
2m
Q(r) = V (r) ψ(r)
~2
D’où,

m +∞ ejk |r−r0 |
Z
ψp (r) = − V (r0 ) ψ(r0 ) dr0 (262)
2π ~2 −∞ |r − r0 |
La solution globale (solution homogène ψh (r) + solution particulière ψp (r)) de l’équation de schrödinger
est :

m +∞ ejk |r−r0 |
Z
ψ(r) = ψh (r) − V (r0 ) ψ(r0 ) dr0 (263)
2π ~2 −∞ |r − r0 |

13. Afin de ne pas alourdir les écritures mathématiques, le signe vecteur est délibérément omis. Ainsi, le vecteur r0 décrit tout
les points où règne le potentiel, d’un autre côté, le vecteur r décrit les points à partir desquels nous observons la fonction d’onde.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 101

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

La solution de l’équation homogène :
h i
∇2 + k 2 ψh (r) = 0 (264)

vaut ψh (r) = ejk r , onde plane d’amplitude unitaire. Ainsi, la solution globale (263) devient :

m +∞ ejk |r−r0 |
Z
jk r
ψ(r) = e − V (r0 ) ψ(r0 ) dr0 (265)
2π ~2 −∞ |r − r0 |
L’équation (265) est la forme intégrale de l’équation de schrödinger. Elle est totalement équivalente à la
forme différentielle, plus familière. Le ”paradoxe” de cette équation intégrale réside dans le fait que nous devons
d’abord connaitre la solution ψ(r) afin de résoudre cette intégrale, or c’est ce que nous cherchons !. Avant de
répondre à cette question, supposons que r0 ≈ 0 cela implique que V (r0 ) ≈ 0 pour des régions loin du centre
de diffusion du potentiel. Autrement dit, pour de longues distances du centre de diffusion nous avons |r| > |r0 |,
ce qui conduit à l’approximation suivante :

ejk |r−r0 | ejk r −jk r0

≈ e (266)
|r − r0 | r
Par voie de conséquence, l’équation (265) devient :

m ejk r +∞
Z
ψ(r) = e jk r
− e−j k r0 V (r0 ) ψ(r0 ) dr0 (267)
2π ~2 r −∞

Il faut bien rappeler que le phénomène physique rattaché à cette équation, est celui d’une onde plane
incidente, ψ(z) = ejki r voyageant dans la direction ~ui , qui rencontre un potentiel de diffusion V (r0 ), produisant
une onde sphérique sortante ejks r0 dans la direction ~us . Comme une image, nous pouvons imaginer la situation
d’une vague d’eau (donc une onde) rencontrant une roche qui va la dévier selon une direction donnée. Le nombre
d’onde k est liée à l’énergie de la particule à travers la relation :
√
Cours complet est disponible sur mon site web : [Link]

2mE
k≡
~
Tenant compte des directions d’incidence et d’émergence, la solution (267) devient :

m ejki r +∞
Z
ψ(r) = ejki r − e−j ks r0 V (r0 ) ψ(r0 ) dr0 (268)
2π ~2 r −∞

Afin de résoudre cette équation intégrale, nous devons passer par l’approximation de Born 14 . Cela consiste
à écrire :

m ejki r +∞
Z 0 0 0 0
ψ(r0 ) = ejki r0 − e−j ks r0 V (r0 ) ψ(r0 ) dr0 (269)
2π ~2 r0 −∞

Substituons ensuite l’équation (269) dans (268), il vient :

m ejki r +∞ m ejki r +∞
Z Z 0 0 0 0
ψ(r) = ejki r − e−j ks r0 V (r0 ) ejki r0 dr0 − e−j ks r0 V (r0 ) ψ(r0 ) dr0 + · · ·
2π ~2 r −∞ 2π ~2 r0 −∞
| {z }
premier approx. de Born

Si l’on tronque cette série à la première approximation de Born, la solution approximative prend la forme
suivante :

m ejki r +∞
Z
ψ(r) = ejki r − e−j ks r0 V (r0 ) ejki r0 dr0 (270)
2π ~2 r −∞

14. Approximation ayant été introduit dans les années 1920 par Max Born, physicien théoricien Allemand.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 102

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Ou encore :

ejki r m +∞
Z
ψ(r) = ejki r − ej (ki −ks ) r0 V (r0 ) dr0 (271)
r 2π ~2 −∞
| {z }
f (θ,φ)

Avec,
m
Z +∞
f (θ, φ) = − ej κ r0 V (r0 ) dr0 , κ ≡ ki − ks et r > r0 (272)
2π ~2 −∞

Désormais la connaissance de V (r) permet la détermination de la solution ψ(r) par le biais de l’équation
intégrale (271). La dépendance angulaire de l’amplitude de l’onde émergente f est portée par le nombre d’onde
κ. Tout la problématique consiste donc à déterminer l’amplitude de l’onde émergente f (θ, φ), ce facteur indique
la probabilité de l’émergence dans une direction donnée θ.

jks
jki

ks

ki

Figure 9: Onde plane incidente (~ki ), diffusée par le centre de diffusion V (r) ayant une porté de l’ordre del.
Cours complet est disponible sur mon site web : [Link]

Le résultat est une onde sphérique sortante (~ks ) ayant une amplitude f (θ, φ) dans l’angle solide Ω.

D’après le schéma, nous pouvons écrire :

κ
2
⇒ κ = ks − ki = 2k sin(θ/2)
sin(θ/2) = (273)
k
Où ki et ks ont strictement la même amplitude que le nombre d’onde k, mais le premier pointe dans la
direction du faisceau incident, tandis que le second pointe vers le détecteur. Cette égalité des amplitudes
découle du caractère élastique de l’interaction entre l’onde incidente et le centre de diffusion V (r). Autrement
dit, cette interaction ne change pas l’énergie de l’onde incidente, d’oú |ki | = |ks | = |k|. L’amplitude de l’onde
émergente f (θ, φ), se réécrit selon :
m
Z +∞
f (θ, φ) = − e−j 2 k sin(θ/2) r0 V (r0 ) dr0 , r > r0 (274)
2π ~2 −∞

Pour une diffusion de faible énergie, autrement dit, si l’onde incidente ki n’est pas trop altérée par le potentiel
V (r0 ), nous pouvons écrire :
2π
E= ~c
λ
|{z}
k

E→0 ⇒ λ→∞ ⇒ k→0 ⇒ e−j 2 k sin(θ/2) r0 → 1

m
Z +∞
⇒ f (θ) ≈ − V (r0 ) dr0 , r > r0 (275)
2π ~2 −∞

Chapitre3 Équations aux dérivées partielles et fonctions de Green 103

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

En considérant un potentiel sphérique et quelque soit l’énergie mobilisée (donc pas forcément dans le cas
des faibles énergies), nous exprimons r0 en coordonnées sphériques. A partir de (272), dans l’argument de
l’exponentielle nous avons un produit scalaire entre deux vecteurs :

κ r0 = ||κ|| ||r0 || cos(θ0 )

Afin de ne pas alourdir les écritures, nous conserverons la notation κ r0 pour les deux modules. Alors :
m
Z +∞
f (θ, φ) = − ej κ r0 cos(θ0 ) V (r0 ) r02 dr0 sin(θ0 ) dθ0 dφ0 (276)
2π ~2 −∞

m
Z 2π Z ∞ Z π 
j κ r0 cos(θ0 )
⇒ f (θ, φ) = − dφ0 V (r0 ) r02 dr0 e sin(θ0 ) dθ0 (277)
2π ~2 0 0 0
" " #π #
m ∞ ej κ r0 cos(θ0 )
Z
⇒ f (θ, φ) = − 2π V (r0 ) r02 dr0 − (278)
2π ~2 0 j k r0 0

m
Z ∞ 
2 sin(κ r0 )

⇒ f (θ, φ) = − V (r0 ) r02 dr0 (279)
~2 0 κ r0
2m ∞
Z
⇒ f (θ) = −
r0 V (r0 ) sin(κ r0 ) dr0 (280)
κ ~2 0
Il suffit maintenant de calculer l’intégrale restante sur la variable radiale r0 . Si nous utilisons un potentiel
de Coulomb simple où V (r) ∝ 1/r, alors nous nous heurtons à un problème car l’intégrale ne converge pas
comme r → ∞. Pour cette raison, un autre potentiel radialement symétrique est introduit, qui est donné par :

e−r/l
V (r) = (281)
r
Où l > 0 est une constante décrivant la portée de l’interaction. Ce type de potentiel est connu sous le
Cours complet est disponible sur mon site web : [Link]

nom de potentiel de Coulomb écranté 15 . Le paramètre l détermine la gamme sur laquelle le potentiel est
influent. Il nous permet d’évaluer analytiquement l’amplitude de diffusion. Ce potentiel est susceptible de
décrire l’interaction électromagnétique entre deux particules chargées s’il y a un écrantage de la force lorsque
la distance entre les particules est plus grande que l’ordre de la distance caractéristique ∝ l. Il est possible
alors d’observer le comportement de f (θ, φ)2 pour un potentiel de Coulomb en laissant l tendre vers zéro.
L’amplitude de diffusion devient :
2m
Z ∞
⇒ f (θ) = − 2 e−r0 /l sin(κ r0 ) dr0 , κ = 2 k sin(θ/2) (282)
κ~ 0

2m κ
 
⇒ f (θ) = − 2 , κ = 2 k sin(θ/2) (283)
κ~ l + κ2
2

2m 2 k sin(θ/2)
 
⇒ f (θ) = − (284)
(2 k sin(θ/2)) ~2 l + (2 k sin(θ/2))2
2

2m 1
 
⇒ f (θ) = − 2 (285)
~ l + (2 k sin(θ/2))2
2

Ainsi, lorsque le paramètre l se rapproche de zéro, on obtient :

2m 1
 
⇒ f (θ) ≈ − (286)
~2 (2 k sin(θ/2))2

15. Ce potentiel s’obtient en superposant une charge ponctuelle positive à l’origine et une charge négative délocalisée au
voisinage de l’origine sur une distance de l’ordre de l.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 104

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Par conséquent, l’intensité de l’onde émergente (ou dispersée) est :
1
I = |f (θ)|2 ∝ 4 (287)
sin (θ/2)
Dans cette analyse, nous avons considéré le cas d’un potentiel d’interaction central V (r) qui dépend seule-
ment de la distance radiale par rapport au centre du potentiel. Il est intéressant de noter que l’intensité
d’interaction |f (θ)|2 est indépendante des variables caractérisant l’onde incidente (intensité, densité, ... etc.)
et le potentiel en question, elle dépend uniquement de l’angle de déviation.
Cours complet est disponible sur mon site web : [Link]

Chapitre3 Équations aux dérivées partielles et fonctions de Green 105

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

IV. Annexe : Compléments sur les orbitales des atomes réels
Nous entendons par l’adjectif réels les atomes ayant plus d’un électrons, ou simplement les atomes poly-
électroniques. En effet, bien que les orbitales atomiques des ions hydrogénoı̈des soient des solutions exactes de
l’équation Schrödinger donc répondant de manière rigoureuse au formalisme de la mécanique quantique, elle
sont néanmoins peu efficace pour d’écrire les propriétés des atomes réels et par extension celles des molécules.

A. Orbitales de type Slater

A cet égard, John Clarke Slater a proposé de garder la même forme générale que celle des orbitales atomiques
hydrogénoı̈des mais toutefois avec un terme radial dépendant uniquement cette fois-ci du nombre quantique
principal. Dans cette démarche, John C. Slater propose également une constante d’écran afin de prendre en
compte l’effet de la distribution de charge des électrons du cœurs. L’expression mathématique des orbitales
atomiques de Slater (ces orbitales sont regroupées sous l’acronyme STO, pour Slater Type Orbital) est donnée
par :

χST O m
n,l,m (r, θ, φ) = Rn (r) × Yl (θ, φ) (288)
Où la partie radiale est :
2n+1
(2 ζ) 2
Rn (r) = (1/2)
rn−1 e−ζ r (289)
(2 n!)
Les mêmes harmoniques sphériques Ylm (θ, φ) des atomes hydrogénoı̈des sont utilisées pour décrire la partie
angulaire. L’exposant ζ (qui se lit zêta en alphabet latin) est un paramètre ajustables, lié à la charge effective du
noyau. La charge nucléaire étant partiellement ”écrantée” par les électrons des couches internes. Cette charge
effective est estimée par les règles de Slater que détaillerons ci-dessous. En outre, le paramètre ζ controle
la largeur de l’orbitale, une valeur élevée de ζ produit une orbitale plus fine et petite valeur produira au
contraire une orbitale plus large. A titre d’exemple pour représenter l’orbitale atomique |1si, nous utiliserons
l’expression :
Cours complet est disponible sur mon site web : [Link]

" #1/2
ζ3
χST O
1,0,0 (r, θ, φ) = e−ζ r × Y00 (θ, φ) (290)
π
Par ailleurs, il est possible utiliser plus d’une STO pour représenter une orbite atomique, comme le montre
l’équation :
h i
−ζ1 r
χST O
2,0,0 (r, θ, φ) = c1 r e + c2 r e−ζ2 r × Ylm (θ, φ) (291)
Dans cette équation, les paramètres ζ1 et ζ2 sont ajustés par une procédure de fitting (les moindres carrés par
exemple). Les coefficients c1 et c2 sont déterminés par un calcul variationnel linéaire qui minimise l’énergie du
système quantique étudié. Dans cette configuration, la fonction ayant la plus grande valeur de ζ tient compte
de la charge près du noyau. Tandis que la fonction ayant la petite valeur de ζ tient compte de la distribution
de la charge à des valeurs plus importantes de la distance du noyau. Cette base de fonctions d’onde est appelée
ensemble, base 16 à double zêta. Nous exigeons que ces fonctions de base couvrent tout l’espace de distribution
des électrons 17 , ce qui signifie qu’elle doivent former un ensemble complet et doivent décrire la même chose.
Par exemple, les harmoniques sphériques ne peuvent pas être utilisées pour décrire une fonction radiale d’un

16. Une base de fonctions d’onde (basis set en anglais) en mécanique quantique est un ensemble
X de fonctions (appelées fonctions
de base) qui sont linéairement combinées pour créer des orbitales moléculaires : ψi = cij ϕj . Par commodité, ces fonctions
j
sont généralement des orbitales atomiques centrées sur les atomes.
17. Cette propriété des fonctions dans l’espace est tout comme la propriété correspondante des vecteurs. Les vecteurs unitaires
(~e1 , ~e2 , ~e3 ) décrivent des points dans l’espace et forment un ensemble complet puisque toute position dans l’espace peut être
spécifiée par une combinaison linéaire de ces trois vecteurs unitaires. Ces vecteurs unitaires sont également appelés vecteurs de
base.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 106

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

atome d’hydrogène car elles n’impliquent pas la variable radiale, mais elles peuvent être utilisées pour décrire
les propriétés angulaires de n’importe quel système dans l’espace tridimensionnel.

Par ailleurs, nous pouvons démonter que l’énergie d’un électron occupant une orbitale de Slater à exactement
la meme forme que celle obtenu en résolvant l’équation de Schrödinger. En utilisant le Laplacien en coordonnées
sphériques, l’équation de Schrödinger devient :
" #
~2 1 ∂ ∂ 1 ∂ ∂ 1 ∂2
   
− 2
r2 + 2 sin θ + 2 2 χST O + V (r) χST O = E χST O (292)
2m r ∂r ∂r r sin θ ∂θ ∂θ r sin θ ∂φ2
| {z }
∇2

En adoptant le système des unités atomiques nous avons : m = ~ = e = 4 π 0 = 1 et afin de simplifier

encore cette expression, posons aussi :

1 ∂ ∂ 1 ∂2
 
sin θ + 2 =Λ (293)
sin θ ∂θ ∂θ r2 sin θ ∂φ2
Nous obtenons alors :
1 1 ∂ ∂ Λ
   
− 2
r2 + 2 χST O + V (r) χST O = E χST O (294)
2 r ∂r ∂r r | {z }
Z∗
−
r
Les orbitales de Slater sont données par :

χST O m
n,l,m (r, θ, φ) = Rn (r) × Yl (θ, φ) (295)
Où la partie radiale dépendant uniquement du nombre quantique n est :
Cours complet est disponible sur mon site web : [Link]

2n+1
(2 ζ) 2 Z∗
Rn (r) = rn−1 e−ζ r avec ξ = (296)
(2 n!)(1/2) n
Il est bien connu que les harmoniques sphériques Ylm (θ, φ) sont des fonctions propres de l’opérateur Λ :

λ Ylm (θ, φ) = −l(l + 1) Ylm (θ, φ) (297)

En substituant (295) dans (294) il vient :

−l(l+1) Ylm (θ,φ)

 
{ z }|
Rn (r) Λ × Ylm (θ, φ)  ∗
 m 
1  Yl (θ, φ) ∂ 2 ∂ Rn (r)  − Z Rn (r)×Y m (θ, φ) = E Rn (r)×Y m (θ, φ) (298)
 
−  r + l l
2  r2 ∂r ∂r r2 
 r

Calculons d’abord la dérivée :

2n+1
∂ ∂ Rn (r) (2 ζ)
  2 h i
r2 = e−ξ r (n − 1)2 rn − 2 n ξ rn + ξ 2 rn+1 (299)
∂r ∂r (2 n!)(1/2)
2n+1
∂ ∂ Rn (r) (2 ζ)
  2 h i
⇒ r2 = rn−1 e−ξ r (n − 1)2 r − 2 n ξ r + ξ 2 r2 (300)
∂r ∂r (2 n!)(1/2)
| {z }
Rn (r)

Chapitre3 Équations aux dérivées partielles et fonctions de Green 107

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

En substituant (300) dans (298) il vient :
" #
1 χST O h i l(l + 1)
ST O Z ∗ ST O
− (n − 1)2
r − 2 n ξ r + ξ 2 2
r − χ − χ = E χST O (301)
2 r2 r2 r
" #
1 −l(l + 1) (n − 1)2 − 2 n ξ − 2 Z ∗
⇒ − + + ξ 2 χST O = E χST O (302)
2 r2 r
" #
1 −l(l + 1) (n − 1)2 − 2 n ξ − 2 Z ∗
⇒ En = − + + ξ2 (303)
2 r2 r
" #
1 −l(l + 1) (n − 1)2 − 2 n ξ − 2 n ξ
⇒ En = − + + ξ2 (304)
2 r2 r
Pour une distance suffisamment loin du noyau, c’est-à-dire r grand, nous pouvons considérer les approxi-
mations :

−l(l + 1) (n − 1)2 − 4 n ξ
−→ 0 et −→ 0 (305)
r2 r
1 1 Z∗ 2
 
⇒ En ' − ξ 2 = − (306)
2 2 n
Le traitement que nous venons d’effectuer nous informe que les orbitales de Slater décrivent ”bien” le
comportement de l’électron tant que celui-ci ”se tient” loin du noyau. Pour les atomes hydrogénoı̈des le
paramètre Z exprime la charge du noyau. Pour des atomes non-hydrogénoı̈des, ce paramètre est la charge
effective Z ∗ ressentie par l’électron. Selon le modèle de Slater, la charge effective ressentie par l’électron é(i)
est déterminée selon :

Zi∗ = Z −
X
Nj σ j avec 1 < j ≤ i (307)
Cours complet est disponible sur mon site web : [Link]

j≤i

Le paramètre Nj représente le nombre d’électrons des couches nj ≤ ni . Dans ce modèle de Slater, l’inter-
action Coloumbienne entre un électron (i) et son noyau est supposée être perturbée ou ”écrantée” par les
électrons (j) des couches intermédiaires. Cet effet ”d’écrantage” de l’intération Coloumbienne électron-noyau
est interprétée en terme d’une constante dite constante d’écran portant le symbole σj . Cette constante dépend
de l’emplacement des électrons j (N − 1 électrons restant) situés entre l’électron i et le noyau. Par conséquent
l’électron i ressentira une charge réduite +Z ∗ |e| au lieu de la charge réelle du noyau +Z |e|.

B. Orbitales Gaussienne
Lorsque des calculs quantiques sur des moléculaires sont menés, il est courant d’utiliser une base composée
d’un nombre fini d’orbitales atomiques, centrées sur chaque noyau atomique de la molécule. Ces orbitales
atomiques sont bien décrites avec les orbitales de type Slater, car les STO décroissent exponentiellement avec
la distance par rapport au noyau, décrivant ainsi précisément le chevauchement à longue distance entre les
atomes. Néanmoins, il a été démontré 18 que les intégrales impliquant des fonctions gaussiennes sont plus
rapides à calculer que celles impliquant des exponentielles de type STO, ce qui se traduit par un gain en terme
de temps de calcul. D’une façon générale, les orbitales gaussiennes s’écriront sous la forme :
l 2(n−l−1) −ζ r 2
χGT O
n,l,m (r, θ, φ) = Nn,l (ζ) r r e × Ylm (θ, φ) (308)
Où Nn,l (ζ) est une constante de normalisation. A titre d’exemple pour représenter l’orbitale atomique |1si,
nous utiliserons l’expression :
3/2
2ζ

2
χGT O
1,0,0 (r, θ, φ) = e−ζ r × Y00 (θ, φ) (309)
π
18. Dans les années 1950, par Frank Boys de l’université de Cambridge au Royaume-Uni.

Chapitre3 Équations aux dérivées partielles et fonctions de Green 108

 M. Samir KENOUCHE

MASTER PHYSIQUE - Toutes options confondues

Notons que pour toutes les foncions de base, seule la partie radiale de l’orbitale change. Les fonctions
harmoniques sphériques sont utilisées décrire la partie angulaire de l’orbitale. Malheureusement, les fonctions
gaussiennes ne reflètent pas précisément la forme d’une orbitale atomique. En particulier, elles sont plutôt
plates que raides près du noyau atomique r → 0. En outre, ces orbitales gaussiennes diminuent plus rapidement
pour les valeurs élevées r → ∞. Afin de remédier à ces limitations, chaque orbitale STO est remplacée par
un certain nombre de fonctions gaussiennes avec des valeurs différentes pour le paramètre ζ. Ces fonctions
gaussiennes forment un ensemble de bases gaussiennes primitives. Des combinaisons linéaires des Gaussiennes
sont formées pour se rapprocher de la partie radiale d’une orbitale STO. Ces nouvelles fonctions sont définies
alors comme desfonctions gaussiennes contractées (CGTO) :
N
2
cζj rl r2(n−l−1) e−ζj r × Ylm (θ, φ)
X
CGT O
χ (r, θ, φ) = (310)
j

Où N est la taille de la contraction et cζj sont les coefficients de la contraction. Les fonctions gaussiennes
contractées les plus simples sont les bases STO-nG. Cette base tente d’approcher ou de mimer les orbitales
STO par la somme de N Gaussiennes.
Cours complet est disponible sur mon site web : [Link]

Chapitre3 Équations aux dérivées partielles et fonctions de Green 109