W g a m s I n s titu u t v o o r á i I k s

Flandern Marino Instituts

'h y sicalia M ag. 21 (1 9 9 9 ) 1 pp 31-56

L’EQUATION DE LÀ CHALEUR ET L’ENSEIGNEMENT DES

METHODES DE SIMULATION NUMERIQUE DE LA
PHYSIQUE
Eric D eleersnijder

In s titu t d ’A s tr o n o m i e e t d e G é o p h y s iq u e G. L e m a îtr e
U n iv e rs ité C a th o liq u e d e L o u v a in , 2 C h e m in d u c y c lo tro n , 1 3 4 8 L o u v a in - la -N e u v e

A lth o u g h [...] m o r e e n g in e e r s now u se " o ff-th e -s h e lf' s o ftw a r e p a c k a g e s , it is s till

i m p o r t a n t to u n d e r s t a n d w h a t is in th e m a n d how to u se th e m p r o p e r ly . It is th e
a u t h o r 's e x p e r ie n c e t h a t o n e c a n o b ta in th is k in d o f u n d e r s ta n d in g o n ly b y a c tu a l ly
w ritin g a t l e a s t s o m e s im p le p r o g r a m s a n d "p la y in g " w ith th e m e th o d s ; o n e d o e s
n o t r e a lly u n d e r s t a n d a m e th o d u n til o n e h a s p r o g r a m m e d it. I h a v e m e t m a n y
p e o p le w ho u s e c o m m e r c ia l p a c k a g e s a n d c o m p la in a b o u t th e ir e x p e r ie n c e with
th e m , i t o f t e n t u r n s o u t t h a t t h e y a r e e ith e r a s k in g a p r o g r a m to d o s o m e t h i n g i t
w a s n o t i n t e n d e d to d o o r t h a t t h e y a re tr y in g to s o lv e a p r o b le m t h a t d o e s n o t
h a v e a so lu tio n (a n ill-p o se d p ro b le m ). A little b a c k g r o u n d in th e fu n d a m e n t a l s o f
how th e m e t h o d w o rks, w h a t it is g o o d for, a n d h o w it s h o u ld b e u s e d w o u ld h a v e
h e lp e d th e s e p e o p le e n o r m o u s ly .
Joel H. Ferziger (1998) (à propos des méthodes numériques de l ’ingénierie)

Sum m ary. N um erically solving partial differential problem s to sim ulate the beh av io u r o f
natural or artificial system s is n o w p art o f the scientific m ethod — along w ith the experim ental
and theoretical approaches — m ainly because the pow er o f the com puters h a s enom ously
increased over the last decades. T his is the reason w hy it is argued that p h y sics students should
be given an initiation to m ethods fo r solving partial differential problem s — preferably during
their third y ear at th e university. It is believed th at basic properties o f fm ite-difference numerical
schem es — i.e. consistency, stability, convergence, conservativity, accuracy, efficiency — may
be introduced and illustrated w hile building m ethods fo r num erically solving sim ple partial
differential equatio n s supposedly know n to all physicists, such as the heat equation.
A ccordingly, several heat conduction problem s are set up in such a w ay that the analytical
solutions o f the continuous problem s and their discrete counterparts are available, rendering it
easy to understand th e properties o f the corresponding solution m ethods. A particular em phasis
is p u t on one-dim ensional p ro b lem s, b u t a few tw o-dim ensional m ethods are also considered.
F inally, the im portance o f advective heat transport processes is m entioned, as w ell as the
specific difficulties relevant to the num erical treatm ent o f these phenom ena. H o p efu lly , these
exam ples w ill be useful to those w illing to introduce physics students to the fundam entals of
num erical sim ulation in physics.

31
I n t r o d u c tio n

D epuis le X V IIèm e siècle environ, on a construit des m o d èles m athém atiques1, c'est-à-d ire des
relations m athém atiques susceptibles de décrire l'état de systèm es naturels ou artificiels, ainsi
q u e l'év o lu tio n de cet état. En d'autres term es, ces sy stè m e s ont été “m is en éq u atio n s” . Ces
dernières so n t souv en t d es équations différentielles o rd in aire s ou aux d érivées partielles. On a
établi de n o m b reu ses propriétés qualitatives utiles des so lu tio n s de ces éq u atio n s, m ais il a
rarem ent été p o ssib le d'en découvrir les solutions analytiques g énérales. D ans la p lupart des
cas, obtenir des so lu tio n s analytiques n'a été possib le que p o u r d es problèm es idéalisés,
découlant d 'h y p o th è se s fortes sur la sym étrie des so lu tio n s, la form e d u dom aine d 'in térêt, la
proxim ité d 'un état d 'éq u ilib re , etc. C es difficultés m athém atiques on t san s doute freiné
l’interaction entre les développem ents théoriques et l'expérim entation ou l'observation.
R echercher une solution approxim ative d'un problèm e m athém atique constitue
fréquem m ent une alternative fructueuse à la déterm ination de la solution exacte. A insi, les
solutions asym ptotiques so n t connues depuis plusieurs siè cle s, m ais on ne p eu t construire de
telles approxim atio n s que dans des conditions bien p récises, qui lim itent drastiquem ent la
construction e t l'em ploi de ces expressions. D 'un autre c ô té , il est quasim ent to u jo u rs possible
de bâtir d es algorithm es num ériques susceptibles de fo u rn ir une approxim ation de la solution
exacte d 'u n pro b lèm e. Cette approche n'a cependant p ris un réel esso r q u 'av ec l'ap p aritio n de
calculateurs su ffisam m en t p uissants (Figure 1), c'est-à-d ire vers la fin des années septante. Les
prem iers o rdinateurs, qui datent des années quarante, étaient capables d 'e ffe c tu e r quelques
opérations arithm étiques p ar seconde. A ujourd'hui, u n e station de travail com m une peut
réaliser, chaq u e seco n d e, de l'o rd re de 1 0 0 m illions de m ultiplications de n om bres com portant
environ 16 ch iffres significatifs. U ne telle puissance de calcul perm et de p ro d u ire des
approxim ations de la solution de problèm es m athém atiques relativem ent com plexes de sorte
qu'on p eut m aintenan t sim uler num ériquem ent avec un h a u t niveau de réalism e le com portem ent
de nom breux systèm es naturels ou artificiels. On peut affirm er que la sim ulation num érique est
m aintenant devenue un acteur à part entière de la dém arche scientifique2, au m êm e titre q u e les
développem ents théoriques et l'expérim entation en laboratoire, ou encore l'o b serv atio n su r le
terrain.
A u jo u rd ’hui, la sim ulation num érique tend à se substituer aux expériences sur des systèm es
en grandeur réelle ou réduite dès que ces dernières so n t p lu s chères q u e la sim ulation — à
con d itio n , bien en ten d u , que la simulation du systèm e considéré soit possible. A insi, pai'
exem ple, on rem place d e plus en p lus les essais en soufflerie p ar des program m es inform atiques
visan t à réso u d re les équations de l ’aérodynam ique — q u e l’on qualifie p arfo is de “ soufflerie
n u m ériq u e”. D ’autre part, il existe des systèm es inaccessibles à l’expérim entation au sens
classique d u term e. L e systèm e clim atique de la T en'e appartient clairem ent à cette catégorie: on

1 L e s P r in c ip ia d e N e w to n , p u b lié s e n 1687, fu ren t san s doute l e p re m ie r grand succès d e la m odélisation

m a th é m a tiq u e , e n p ro u v a n t c la ire m en t la capacité des m odèles m a th é m atiq u es à rep résenter e t, d ans une certaine
m e su re, e x p liq u e r d e s p h é n o m è n es naturels.
2 A p ro p o s d e la s im u la tio n n u m é riq u e , on ap p ré ciera c es m o ts de R o b in so n (1987): “ T h e use o f su ch s im u la tio n
in s c ie n tific research (n u m erical e x p erim e n tatio n , sen sitiv ity and p ro c e ss studies, e tc .) is th o u g h t b y m a n y to
re p re sen t th e firs t m a jo r ste p fo rw ard in th e b asic scien tific m e th o d sin c e th e sen venteenth c en tu ry . S cien ce is
no w a trip a rtite e n d e a v o r w ith sim u la tio n a d d ed to th e tw o classical c o m p o n e n ts. E x p e rim en t an d T h e o ry ” .

32
ne p eut “m ettre un tel systèm e en éprouvette” pour étudier e t com prendre son fonctionnem ent,
m ais on p eut sim uler son com portem ent à l’aide de m oyens inform atiques et p ratiq u er alors des
“expériences num ériq u es” .

puissance des calculateurs

1940 1960 1980 2000

t: année
F ig u r e 1. C o u rb e rep résentant l ’évolution d epuis la fin d es années q u arante de la
vitesse approxim ative du p ro cesseu r le p lus puissan t de l ’époque, exprim ée en
nom bre d ’opérations arithm étiques p ar seconde (selon H ack, 1992).

N om bre de sim ulations num ériques rep o sen t su r la production d 'ap p ro x im atio n s discrètes
de la solution de problèm es différentiels. Les algorithm es correspondants so n t b asés su r des
schém as num ériques d ont la m ise au point n 'est généralem ent pas triviale. S 'il existe
actuellem ent d es logiciels com m erciaux — com m e Mathematica ou M atlab — perm ettant de
fournir des solutions num ériques d e problèm es différentiels aux d érivées ordinaires — san s que
leur utilisateur doiv e s'im p liq u er dans les détails des schém as num ériques — , la résolution
num érique de problèm es différentiels aux dérivées partielles exige souvent la m ise au point d 'u n
algorithm e spécialem ent conçu p o u r le problèm e considéré. Le scientifique ne p eu t donc s 'e n
rem ettre à un logiciel qui prendrait en charge les “détails techniques” de la résolution num érique
du problèm e. E t m êm e qu an d u n tel logiciel existe, l'u tiliser valablem ent dem ande des
com pétences num ériques significatives. En d ’autres term es, pratiquer la sim ulation num érique
— et, en particulier, la résolution num érique de problèm es différentiels aux d ériv ées partielles
— ne s ’im p ro v ise pas.
P our sim uler n u m ériq u em en t u n systèm e, il faut d'abord établir les équations q u i régissent
l'état de ce systèm e. E nsuite, il faut bâtir une ou plusieurs approxim ations discrètes de ces
équations et, finalem ent, im plém enter sur un ordinateur les program m es co rresp o n d an ts. On
peut pen ser que le program m e de la licence en physique des universités fran co p h o n es belges

33
couvre actuellem ent très bien la prem ière étape, à savoir la m ise en équations du sy stè m e , ainsi
que l'établissem ent des propriétés qualitatives des so lu tio n s, la découverte d e solutions
analytiques d an s des cas idéalisés et, dans une m oindre m esu re, la construction
d'app ro x im atio n s asym ptotiques. P ar contre, il est regrettable que, d an s la p lupart d es cas, les
étudiants ne connaissen t quasim ent rien des m éthodes num ériques visant à o btenir des
approxim ations discrètes des problèm es différentiels — su rto u t aux dérivées p artielles. Il serait
peut-être opportun de com bler cette lacune, afin que n o s licenciés en phy siq u e so ien t mieux
arm és p o u r obtenir des em plois dans les dom aines — toujours p lus nom breux — de la
recherche fondam entale ou appliquée où la simulation n u m ériq u e jo u e un rôle im portant.
Initier les étudiants aux m éthodes de simulation n u m érique de la p h y siq u e — et, en
particulier, à la résolution num érique d'équations différentielles aux dérivées partielles —
n'im plique p as nécessairem ent une lourde charge d 'en seig n em en t. En e ffet, en prem ière
licence3, on p o u rra it se fixer l'objectif d'introduire et illu strer la plupart d es notions de base de
la technique d es différences finies 4 en traitant quelques p ro b lèm es différentiels élém entaires qui
concernent tous les ph y sicien s. C i-dessous, on va m o n tre r que la résolution de p roblèm es de
conduction de chaleur peut aisém ent aider à atteindre ce but.

E q u a tio n d e la c h a le u r u n i-d im en sio n n e lle

On considère le problèm e classique de conduction de ch a le u r dans une baiae isolée — sau f en

ses extrém ités. La longueur de la barre v au t 2L, et la p o sitio n de chaque point de cette barre est
repérée p a r la coord o n n ée x , telle que - L < x < L . O n divise la barre en 2 K segm ents de
longueur ég ale A x = (2 L )/(2 K ) = L / K . On identifie ch aq u e segm ent à l'a id e de l'in d ic e k , qui
p eut prendre les valeurs suivantes (Figure 2):
k = - K + l , -K + 2 ,... K - l , K . (1)

Le centre d u segm en t associé à l'indice k est situé à la coordonnée

xk = (k -\/2 )A x . (2)

L es extrém ités d u segm ent k sont donc situées en x k+l/2 = (k -1 /2 ± 1/2) A x . La tem pérature
évo lue d ans le tem p s, que l'on note t. On identifie des instants su ccessifs, séparés pai' la durée

3 L a p re m iè re lic e n c e p o u rrait ê tre le m o m e n t idéal p o u r in itie r les é tu d ia n ts à la ré so lu tio n nu m é riq u e d ’équations

différentielles aux d érivées p a rtie lle s. E n candidature, la m a tu rité scientifique — n o ta m m e n t les bases
m ath ém atiq u es e t la cap a c ité à é ta b lir des re latio n s fructueuses e n tre les m atières de différents co u rs — est
in su ffisa n te . E n sec o n d e lic e n c e , l ’in itiatio n au x m éthodes dont il e s t question ici v ien d rait tro p tard p o u r q u ’un
m é m o ire y fa is a n t la rg e m e n t a p p el p u isse ê tre en trep ris. Et, c o m m e le m ém o ire c o n d itio n n e so u v en t le c h o ix des
élu d es et re c h erch e s d e tro isiè m e cy cle — o u, parfois, l ’obtention d ’ u n e m p lo i dans l’in d u strie — , ne p a s ê tre en
m esu re d e p ra tiq u e r v a la b lem e n t la sim ulation num érique en s ec o n d e licence pourrait re v e n ir à y re n o n c er de fa c to
p en d an t u n e p é rio d e sig n ific a tiv e d e la carrière.
4 L e s “ différences fin ies” s o n t u n e fa m ille d e m éthodes n u m ériq u es c ouram m ent e m p lo y é es p o u r résoudre des
p ro b lè m e s d ifféren tiels. S c h é m a tiq u em en t, o n p e u t dire que les a u tre s fam illes sont le s “ é lé m en ts fin is” e t les
“ m é th o d e s s p e c tra le s ” . S i le s fro n tiè re s e n tre ces trois cla sses de m é th o d e s so n t plutôt flo u e s, il e s t p a r c o n tre très
c la ir q u e les d ifféren c e s fin ies co n stitu en t la catég o rie la p lu s s im p le à m e ttre en o e u v re — e t donc à enseigner.
C e p e n d a n t, le s c o n ce p ts d e b a se , c o m m e la p récisio n , stab ilité, la c o n v erg en c e , e tc ., des sc h é m a s n u m é riq u e s sont
valab les p o u r to u te m éth o d e. D o n c, les e n seig n er dans le cadre d e la m éthode des différences fin ies n e p e u t que
faciliter l'ap p ren tissag e u lté rie u r d 'au tres fam illes d e m éthodes num ériques.

34
A t , à l’aide de l'ind ice
/ = 0 , 1, 2 , . . . , (3)
de telle sorte que
t[ = l A t . (4)

P our toute variable a , on utilise les indices k et I com m e suit:

o k = a(tr x k ) . (5)

k-l lk k+ 1
♦
vk - m k+l/2

Ax

F igu re 2. R ep résen tatio n schém atique du Æ-ème seg m en t de barre, avec la position
des tem pératures e t flux de chaleur.

L 'énergie interne com prise dans le ¿-èm e segm ent est définie com m e l’intég rale sur le
segm ent considéré du produit de la m asse volum ique p , de la chaleur spécifique c et de la
tem pérature T. L a barre étant isolée, le prem ier principe de la therm odynam ique indique que
seuls les flux de chaleur aux extrém ités du segm ent considéré sont susceptibles de faire varier
son énergie interne. A insi, si 0 désigne le flux de chaleur transitant d an s la barre — associé au
phénom ène de “cond u ctio n ” — le bilan d'énergie du A-ème segm ent s'écrit
**+1/2 **+1/2 '/+1
J p c T ( t l+ v x ) d x = J p c T ( t r x ) d x + j [<p(t,xk_m ) - <p(t,xk + m )] d t . (6)
S elon la loi de F ou rier, le flux de chaleur tend à transporter de la chaleur de la rég io n la plus
chaude à la région la plus froide e t p eu t être param étrisé com m e suit:
dT
<P = - K - r - , (7)
dx

où k: représente la conductivité therm ique.

Les équations (6 ) et (7) so n t les ingrédients essentiels à l'étab lissem en t de l'équation
différentielle aux dérivées partielles régissant l'évolution de la tem pérature — com m uném ent
appelée “équation de la chaleur” — et d'une approxim ation discrète de cette m êm e équation.
A partir de ce p o in t, p o u r sim plifier les développem ents m athém atiques, on su p p o sera que
la m asse volum ique, la chaleur spécifique et la conductivité therm ique so n t co n stan tes. Par-
conséquent, la diffusivité de tem pérature,

A = JL (8)
pc

35
est égalem ent constante.
E n introduisant les développem ents de Taylor
i dmT \ (* -* * )"
T ( t , x ) = T ( t ,x t ) + £ (9)
m= 1 \ d x ' /

et

c - f,)"
(p(t,x) = <p(tr X) + 2 (10)
»1=1
/n!

d an s ( 6 ), on o btient, ap rès q uelques m anipulations,

Tl,k = Ql.k+l/2 h k - 1/2 + 0 (A , ^ 2 -v

pc (H)
Ar Ar ’
Si l'on augm ente indéfinim ent le nom bre de segm ents, l'ex p ressio n (11) s'a p p liq u e alors “en
chaque point” de la barre. A insi, à la limite A x / L -» 0 o u K (11) conduit à l'équation
continue
dT _ <90
( 12 )
^ dt dx

En intégrant (12) sur le dom aine d'intérêt, on o b tien t la relation rég issan t l’é v o lu tio n de
l'énergie interne de la barre,
+L
E (t) = J pcT dx , (13)
-L

a savoir
0
E (tb) = E (ta) + J [(f)(t,x = - L ) - <p(t,x = + L )] d t , (14)

où tQ et tb sont des instants arbitrairem ent choisis. Ceci m ontre que l'énergie interne de la barre
ne varie q u 'en fonction d u bilan des flux therm iques en ses extrém ités. C 'e st “ physiquem ent”
facile à com prendre: la barre étant isolée, la chaleur ne p e u t en “sortir” ou y “ entrer” que pai- ses
extrém ités! P ar co nséq u en t, si les deux extrémités de la barre sont égalem ent isolées — aucun
flux de ch aleu r ne les traverse — , alors l'énergie interne to tale de la barre reste constante.
S ubstituant la loi de F ourier (7) dans l'équation de bilan de chaleur (1 2 ), utilisant ( 8 ), on
obtient l'équation de la cha leu r uni-dim ensionnelle,

— = A— (15)
dt dx2

qui est le prototype des équations paraboliques5. Si l'on adjoint à cette équation la distribution
initiale de tem pérature et des conditions aux limites adéquates, on obtient un problèm e

5 H irsc h (1 9 8 8 ) (sectio n 3 .2 .3 ) fa it p a rtie des a u teu rs qui ex p o sent les critères pe rm e tta n t de ranger u n e équation
différentielle au x dérivées p a rtie lle s q u asi-linéaire donnée dans la c atégorie ellip tiq u e, parab o liq u e, hyperbolique,
ou d a n s u n e c la sse h y b rid e e n tre d e u x d es trois fo rm es citées p récédem m ent.

36
différentiel aux dérivées partielles d o n t la solution est la distribution de la tem pérature dans la
barre à tout instant postérieur à l'in stan t initial. On peut donc prévoir l'év o lu tio n du champ
therm ique dans la barre — si l'on est en m esure de résoudre le présent problèm e différentiel.

E q u a tio n d is c r è te

L a plu p art d es problèm es différentiels aux dérivées partielles n 'o n t pas de solution analytique
générale connue. P a r contre, on p eu t chercher des approxim ations num ériques de leur solution
p ar la m éthode des différences finies6. Il faut alors renoncer à déterm iner la valeur de la solution
en chaque p o in t e t se contenter d e l'estim er en un nom bre limité de p o in ts. P o u r ce faire, la
relation (11) constitue un point de départ utile. En effet, en négligeant les term es d 'o rd re
supérieur en les incrém ents tem porel e t spatial — 0 ( A t,A x 2) — , on obtient l'algorithm e discret

At 0 /,*+i/2 “
L u - Tu to • <16>

où chaque variable m unie d 'u n astérisque sym bolise l'équivalent discret de la variable continue
correspondante — qui est dépourvue d 'astérisque. L 'ex p ressio n (16) n 'e st évidem m ent pas la
seule expression discrète susceptible de fournir une approxim ation de la solution de ( 1 2 ), mais
elle est sans doute la plus sim ple.
L 'énergie intern e de la solution discrète est
K

E] = X P cTl,I a * ■ ( 1?)
k= -K +1

E n com binant (16) et (17), on ob tien t l'équivalent discret de la relation (14):

e i = e ; + L C * - 0 A i - c «
b
1=1

où la e t lb sont des nom bres entiers d o n t la valeur peut être choisie arbitrairem ent. D onc, p o u r la
solution discrète, com m e p o u r la solution continue, seuls les flux de chaleur trav ersan t les
extrém ités de la barre so n t susceptibles de m odifier l'énergie interne: on dit q u e le schém a
num érique (16) est conservadf. Ce résulat ne dépend en aucune m anière de la façon d ont on
estim e les flux de chaleur à l'in térieu r de la barre. A insi, par exem ple, cette propriété resterait
valable si on évaluait les flux de chaleur à l'intérieur de la barre à l'aide d'un générateur de
nom bres aléatoires — une option de m odélisation qui ne fournirait évidem m ent pas la meilleure
approxim ation nu m érique de la solution...
En m anipulant le développem ent de T aylor de la tem pérature

> rï
( r)m T 'i 0 - x k±mr
n u x ) = T ( t ,x t ± i a ) + x — ( 19)
m= 1 V ° X J m!
’¿•± 1/2

on p e u t obtenir une approxim ation discrète sim ple de la loi de Fourier,

’ Ici, o n n e fa it p as d e distin ctio n e n tre “ d ifféren ces fin ies” et “ volum es fin is”.

37
 ±T,
l,k± 1 + T,l,k
I ,k ± l/ 2
= — K (2 0 )
Ax
E n introduisant cette dernière dans l'algorithm e (16), il v ie n t
XAt
= T * + ( 2 1 a)
A x 2<-T, M l - 2 T û + T ù - J
l+ lk !,k

ou, une form e équivalente,

T -2 T +T
Tl+\,k Tl,k X I M l l'k l 'k ~ l ( 2 1 b)
At A x2
C et algorithm e constitue une discrétisation élém entaire de l'équation de la chaleur (1 5 ), qui
perm et d 'estim er la tem pérature à l'in stan t i /+1 aux p o in ts x k si la tem pérature est connue à
l'in stan t t¡ aux p o in ts appropriés. A insi, si l'on d isp o se de conditions initiales et aux limites
adéquates, on peut actualiser progressivem ent la tem pérature aux po in ts x k du dom aine
d'intérêt.
En m anipulant les séries de T aylor adéquates, on p e u t ré-écrire (21) sous la form e

- 271*, + 71 dT* ( dLr

i2
T l+ l ,k T i,k l,k +1 l,k ‘W - 1
- X - X (22)
At A x¿ dt dx¿ -Uk
l.k I,k

où e ¡ k désigne l'erreur de troncature, c'est-à-dire la différence, au point d e discrétisation de la

tem pérature (t¡,x k), e n tre l'équation aux différences (2 1 b ) e t l'équation d e la chaleur,
2 T* f d4T*\
e l,k
,. - — i i i A t ------- A x 2 + 0 ( A / 2 ,A r 4 ) . (23)
o
d t 7 I,k 12 d x ¿
l,k

Cette dernière tend vers zéro à m esure que l ’on dim inue le pas de tem ps A t et l ’incrém ent spatial
A x . On dit que le schém a (21) est com patible avec l'éq u atio n de la chaleur (1 5 ), et est du
prem ier ordre de précision selon la dim ension tem porelle et du second ord re selon l'esp ace. Si
l ’erreu r de troncature n 'é tait pas au m oins du prem ier o rd re en chaque incrém ent d es variables
indépendantes, alors cette erreur ne tendrait pas vers zéro à m esure que l'on dim inue ces
incrém ents, indiquant que l'ex p ressio n discrète considérée, représenterait, en fait, une
discrétisation d 'une éq uation continue différente.
E n ex am in an t u n iq u em en t l'errreur de troncature, on pourrait croire qu e, p o u r augm enter la
précision de l'approx im atio n d iscrète d e la solution, c'est-à-dire faire q u 'elle soit p lu s p roche de
la solution exacte, il su ffit d 'au g m en ter l’o rdre de précision du schém a et/ou de dim inuer les
incrém ents d es variab les indépendantes. B ien q u ’elles sem blent en accord avec le bon sens le
p lu s élém entaire, ces sug g estio ns ne sont que partiellem ent correctes: l'erreur de troncature ne
représente p a s une m esu re très fia b le de la précision d u schém a num érique considéré. Deux
raisons com plém entaires expliquent cela:
— L 'erreur d e troncatu re com prend des dérivées de l'in co n n u e dont l'o rd re est le plus souvent
supérieur à celles q u e l'o n trouve dans l'équation continue. Com m e on n 'a généralem ent que
peu d 'inform ations quantitatives ou sim plem ent qualitatives sur ces d ériv ées, la seule
inform ation pertinen te q u e l'on tire de l'analyse de l'erreur de troncature est la puissance la plus

38
petite à laquelle est élevé chaque incrém ent des variables indépendantes — ici, il s'a g it de la
prem ière puissance pour l'incrém ent tem porel et de la seconde pour le pas spatial. O r, u n e telle
inform ation se révèle en définitive peu utile car, dans un processus de résolution n u m ériq u e, les
incrém ents des variables indépendantes ne deviennent ja m a is arbitrairem ent petits7.
— L 'erreur d e troncature est une q u antité à caractère local, car on l'évalue à l'aide des v aleu rs de
la solution prises uniquem ent d an s le voisinage du point considéré de l'esp ace d e s variables
indépendantes. P ar conséquent, cette quantité ne do n n e que peu de renseignem ents, voire
aucun, su r la propagation et l'am plification éventuelle de la différence entre la solution exacte et
son approxim ation num érique.
C 'est parce q u 'elle p eu t être évaluée assez facilem ent et, surtout, sans connaître la solution
ex acte de l'équation continue que l'erreur de troncature est très attrayante. M ais, com m e on vient
de le dém ontrer, cette quantité est peu pertinente et souvent trom peuse. Il faudrait p lu tô t analyser
la convergence du schém a, c'est-à-dire le com portem ent de l'écart entre la solution exacte et son
approxim ation num érique en fonction d es valeurs des incrém ents des variables indépendantes.
P o u r ce faire, il faudrait d isp o ser de la solution exacte. M ais, si celle-ci était c o n n u e, pourquoi
alors en calculer u n e approxim ation discrète? O n touche ici du doigt la difficulté m ajeure
inhérente à tout p rocessus de résolution num érique d'une équation différentielle aux dérivées
partielles: estim er la qualité de la solution num érique est très difficile. Et, puisque l'évaluation de
la solution discrète est m alaisée, déterm iner la façon de l'am éliorer l'est en co re m oins...
L a résolution num érique d'éq u atio n s différentielles aux dérivées partielles n 'est cependant
p as un jeu de hasard. D ifférentes m éthodes existent pour inférer certaines propriétés de
l'approxim ation num érique choisie. P arm i celles-ci, l'an aly se de la stabilité du schém a — poul­
ies équations d'évolution — jo u e un rôle essentiel. Il s'ag it d 'em p êch er que la solution
n u m érique, p o u r des raisons indépendantes de la nature m êm e du p ro cessu s rep résen té par
l'éq uation continue, p résen te des m odes susceptibles de croître indéfinim ent. A insi, un schém a
d 'un ordre de précision élev é au sens de l'erreu r de troncature m ais instable n’a aucune chance
de fo u rn ir une approxim ation valable de la solution exacte, alors q u 'u n schém a d 'u n o rdre de
précision plus faible m ais stable sera très certainem ent p lus utile. P our les problèm es
différentiels aux dérivées partielles vérifiant des conditions m alheureusem ent très restrictiv es, le
théorèm e de Lax (R ichtm eyer and M orton, 1957) stipule que tout schém a com patible et stable
est égalem ent convergent. P our la plupart des problèm es à caractère réaliste, ce théorèm e ne
s'applique pas. On s'y réfère cependant, au m oins im plicitem ent, de façon heuristique: p o u r tout
schém a num érique, on se doit de dém ontrer la com patibilité avec l'équation continue e t analyser,
dans la m esure du possible, la stabilité — ce qui exige souvent une ou plusieurs sim plifications
du schém a, qui com prennent la linéarisation si celui-ci est non-linéaire.

7 L a lim ita tio n de la m é m o ire e t de la p u issa n c e de c alcu l des ordinateurs o b lig e so u v en t le m o d élisateu r à
e m p lo y e r des p a s de te m p s e t d'espace q u i so n t lo in d'être beaucoup p lu s p e tits que les te m p s e t longueurs
c ara c té ristiq u e s d es p h é n o m è n es à sim u ler. A insi, la “ g ro ssièreté” relative du m a illa g e e m p lo y é c o n stitu e l'une des
e x c u se s les p lu s c o u ra m m e n t in v o q u ées p o u r ju s tifie r la p iè tre qualité d'une so lution n um érique.

39
R é so lu tio n d 'u n p ro b lè m e d e c o n d u ctio n de c h a le u r

On va m aintenant illustrer la discussion de la qualité d 'u n e approxim ation num érique en

poursuivant l'étude d u problèm e de conduction de chaleur entam ée p lus haut. P our cela, on doit
prescrire des conditions initiales et aux lim ites. On co n sid ère que les deux extrém ités de la barre
sont isolées. Par conséquent, à tout instant, les conditions aux lim ites

X — = 0 (24)
~ d x î= ± L

s'appliquent. L e profil initial de tem pérature est en “ m arch e d'escalier” :

T ( t = 0 ,x ) = T+ , pour 0 < x < L , (25a)

T ( t = 0 ,* ) = T , pour - L < x < 0 , (25b)

où T + e t T_ sont des constantes. L e problèm e à résoudre se com pose donc de l'éq u atio n (15) et
des conditions (24)-(25).

T a b le a u I. D éfinition d es problèm es adim ensionnels continu et discret de

conduction de chaleur à résoudre.

problèm e continu problèm e discret

dom aine d'intérêt

t > 0 , -1 < * < +1 n = 0 , 1 ,2 ,... ; - K + 1 < k < K

équation d'évolution
rJ~1* 'T'* a 'T ’ * .
dT_ _ d^T_
l+hk ~ l,k l,k+l ~ À 1l,k l l,k - 1
dt ~ dx2 At A x2

condition initiale

T(t=0¿c) = -1 si -1 < x < 0 T 0 k = - 1 si - K + 1 < A: < 0

T (t= 0 jc) = +1 si 0 < a: < +1 Tg k = + 1 si 1 < £ < +K

conditions aux lim ites

T - T
l lt± K +1 l l,± K
= 0
Ax
(S L ■»

Le caractère général de la solution augm ente et les expressions m athém atiques deviennent
p lu s sim ples si l'on rend toutes les variables adim ensionnelles. O n définit ces dernières com m e
suit:

40
 Aí £ T - ( T ++ T _ )/2
(26)
] } 'L (T+ - T _ ) / 2

O n répète la m êm e opération p o u r les variables de l'appro x im atio n discrète de la solution. En

particulier, le pas de tem ps e t le p as spatial adim ensionnels s'écrivent
XAt Ax 1
(27)

A partir de ce p o in t et ju s q u ’à n o u vel ordre, on considérera uniquem ent les variables

adim ensionnelles et, p o u r sim plifier les notations, o n omettra les p rim e s q u i caractérisent ces
variables. O n p eut alors reform uler le problèm e continu et so n pendant discret com m e indiqué
d ans le T ableau I.

0 .5
i= 0.5

1
1 -0 .5 0 0 .5 1
c o o r d o n n é e s p a tia le : x

F ig u r e 3. R eprésentation de la solution exacte (2 8)-(29) d u problèm e de

co nduction de ch ale u r à d ifféren ts instants, à sav o ir t = 0 , 0.0 1 , 0.1 e t 0.5.

L a solution du problèm e continu p o sé dans le Tableau I s'ex p rim e com m e une som m e
in finie de “m odes d e F ourier” (F igure 3):

T ( t ,x ) = ]T a n e v sin (p n x ) , (28)
n= 1

où a n, £n et p n représentent, p o u r le rc-ième m ode, l'am plitude, le co efficien t d 'am o rtissem e n t et

le nom b re d'onde. C es param ètres v alen t
/ V
2
(\ a n ’, e n ,’ “a n /) = f i n - m f n 2 ,{ n -U 2 )n (29)
(rt-l/2 )7 T

41
 En faisant appel à la théorie des équations aux différences (voir, p a r exem ple, B ender and
O rszag, 1978) on obtient 8 la solution analytique du problèm e discret du T ableau I,
K _ *

Tlk = 51 < e ~e" ' • (30)

n= 1

Cette série, au contraire de celle constituant la solution continue, ne com prend q u 'u n nom bre
fini de m odes, K en l'occurence. Si le nom bre d'onde de la solution discrète équivaut à celui du
m ode correspondant d e la solution continue,
f.i* = ( h - 1 / 2 ) k = ß n , pour n < K , (31)

l'am plitude et le coefficient d'am o rtissem en t de la solution discrète ne so n t pas égaux aux
param ètres correspondants de la solution continue. La lo n g u eu r d'onde la plus courte v au t donc
2 A x / ( l - A x / 2 ) — qu i tend vers 2 A x si le n om bre de m ailles tend vers l'infini.

1.5

1.4

d .3
ö

®C1.2

1.1

1 20 40 60 80 100
ordre du mode: n

F ig u re 4. R ap p o rt de l'am p litu d e des m odes, a j a n , p o u r K = 5 (+), K = 10 (o),

K = 50 (x) e t * > 1 0 0 (•)•

L 'am plitude d'un m ode de la solution discrète,

a = ---------, * , -pr:— , (32)

n . (/7-l/2)7T
K s ín --------------
2K

est toutefois asym ptotique à a n si le nom bre de m ailles, K = 1 /A x , est élevé et si l’o rdre relatif,
n /K , du m ode considéré reste petit. En effet, on a

8 S 'il e st rare q u e l'on p u isse tro u v e r la so lu tio n an aly tiq u e d'un p ro b lè m e d ifférentiel a u x d é riv é e s p a rtie lle s, il est
v ra im e n t e x cep tio n n el q u e la so lu tio n an aly tiq u e d 'un problèm e d iscret a ssocié soit d isp o n ib le . C 'e st évidem m ent
p a rc e q u e ces d e u x so lu tio n s a n a ly tiq u e s s o n t c o n n u es — e t re la tiv e m e n t sim p les — q u e l'on a c h o isi de traiter le
p ré sen t p ro b lè m e de c o n d u ctio n d e chaleur.

42
 a * = an + p 0ur n/ K ~ ^ 0 e t K —>oo . (3 3 )
n 12 K
A insi, raffiner la discrétisation spatiale — en augm entant K — am éliore la v aleur d e l'am plitude
d es m o d es d'ordre relativem ent peu élevé. Par contre, q uoi qu 'o n fasse, l'am plitude des m odes
d 'o rd re proche de K reste relativem ent éloignée de celle d es m odes co rresp o n d an ts de la
solution continue (F igure 4).

K = 10 /C = 1 0

1
1

■4
1 5 10 1 5 10
ordre du m ode: n ordre du m ode: n

F ig u r e 5. R eprésentation du facteur d'am plification des m odes n u m ériq u es, y*

(panneau de gauche), et du rapport du facteur d'am plification des m odes
num ériques au facteur d'am plification correspondant de la solution exacte, y* / y ,
p o u r AT=10, d an s le cas où A t / A x 2 - 0 . \ (+), A r/A x 2 = 0 .2 (0 ), A r /A x 2 = 0 .4 (x) et
A t/A x 2= 0 .6 (carré). Les deux prem iers cas co rrespondent à des m odes num ériques
stables, le troisièm e représente des m odes stables m ais d ont certains oscillent, tandis
q u e le dernier est un ex em p le d'instabilité.

P o u r étudier la stabilité des solutions, il est utile d'exam iner le “facteur d ’am plification” , qui
est le facteur par lequel le n-ièm e m ode de la solution est am plifié au co u rs d 'u n p as de tem ps.
A in si, p o u r le /¡-ième m ode de la solution continue et de la solution discrète, le facteur
d'am plification est défini p ar
, *. , -E Al -E * A l.
(ïn ’ 7„) = (e " ,e " ) . (34)

O n v o it sans peine que

0 < y n = e -(n -U 2 f^ A , K { ' (35)

D o n c, aucun m ode de la solution exacte ne croît à m esure que le tem ps p asse. P o u r que le
schém a num érique soit sta b le, ü est indispensable que cette propriété soit vérifiée pai- tous les
m o d es num ériques. Pai' conséquent, on exige que

43
 v i = l - |p r [ l - e o s O Ç A x ) ) (36)

vérifie l'inégalité

M * 1 - (37)

ce qui im plique que le pas de tem ps e t l'incrém ent spatial d o iv en t être tels que

Al < i______ (38)

A. 2 - ! _ c o sí £ d ¿ ! ) £ •
K

Ce résultat est essentiel. Il indique en effet que réduire le pas spatial A x — d ans l'e sp o ir
d'augm enter la précision de l'approxim ation num érique — sans réduire de faço n ad éq u ate le pas
de tem ps p eut très bien conduire à l'instabilité du schém a, c'est-à-d ire détruire totalem ent sa
précision. O n ne p o u va it p a s déduire ceci de la seule analyse de l'erreur de troncature!

1 Af = 0.01
Af = 0.005

0.5 0.5

2 - 0 .5 -0.5
ex
% -1

1 -0.5 0 0.5 1 1 -0.5 0 0.5 1

4
Af = 0.02 At= 0.025

2 - 0 .5 2
■a)
Cl

■4
1 -0.5 0 0.5 1 -0.5 0 0.5 1
c o o r d o n n é e s p a tia le : x c o o rd o n n é e s p a tia le : x

F igu re 6. C om paraison de la tem pérature initiale (tra it pointillé), de la tem pérature

exacte (trait plein) en f = 0 . 2 et de son approxim ation num érique (o) en t= 0 .2 , poul­
i e 5 e t d ifféren tes v aleu rs du pas de tem ps Àt. T o us les schém as so n t stables, sauf
celui correspondant au panneau inférieur droit.

44
 L e facteur d'am plification j n est strictem ent positif p o u r tous les m odes, ce qui im plique
q u 'au cu n m ode n 'o scille dans le tem ps. P o u r qu'il e n soit égalem ent ainsi d e s m odes
n um ériques, il faut leu r im poser une condition plus restrictive q u e la condition de stabilité (38),
à savoir

— < 1 /2 (3 9 1
A* 2 - , _ ' (39)
K
L a F igure 5 illustre le fait que le facteur d'am plification des m odes num ériques d 'o rd re
élévé, c'est-à-d ire proche de K , est, largem ent erroné quelle q u e soit la finesse relative de la
discrétisation spatio-tem porelle.
L e coefficien t d'am ortissem ent des m odes discrets v au t

= ■¿ M 1 - i ? [i ' °os(/i" ^ 4 ■ (4o)

u n e expression qui est com plexe si l'arg u m en t du logarithm e est n égatif. S a u f si le schém a est
in stable, to u s les m odes discrets sont am ortis. P ar conséquent, la différence entre la solution
e x acte e t son approxim ation par un algorithm e stable n 'est jam ais très grande en valeur absolue
p o u r le problèm e qui nous occupe (F igure 6 ). D e toute faço n , on p eu t m ontrer que le schém a
con sidéré ici est convergent dès lors que ses param ètres le ren d en t stable. La dém onstration
form elle de cette propriété so it du cadre de cette note, m ais l'argum entation peut être ébauchée
facilem ent: l'am plitude e t le coefficient d'am ortissem ent des m odes discrets d 'o rd re relativem ent
p eu élevé so n t d'au tan t p lus corrects que les incrém ents spatiaux et tem porels tendent vers zéro
— en p réserv an t la stabilité du schém a — , tandis que l'erreu r su r les m odes d 'o rd re élevé
im porte p eu car ces m odes sont très am ortis.

A u tr e s sc h é m a s n u m é r iq u e s

C 'e st l'ex isten ce d es solutions analytiques continues et discrètes q u i a rendu po ssib le la

d iscu ssion approfondie ci-dessus. L a p lu p art des problèm es différentiels aux dérivées partielles
so n t plus difficiles de telle sorte qu e, généralem ent, aucune solution exacte n 'e st disponible. D
e st cep en d an t souhaitable de disposer d'un outil p o u r déterm iner, au m oins approxim ativem ent,
certaines propriétés des schém as num ériques — en particulier leur stabilité — avant m êm e leur
im plém entation su r ordinateur. C et outil est l'analyse de von N e u m a n n , d o n t on peut
co m p ren d re les bases en se référant à la solution du problèm e résolu ci-dessus. Cette approche
consiste à exam iner le com portem ent de solutions de la form e “m odes de F ourier” d'éq u atio n s
discrètes linéaires à coefficients constants9. Les longueurs d 'o n d e des m odes p euvent prendre
toutes les valeurs supérieures à 2 A x . A insi, en utilisant de nouveau d es variables
d im ensionnelles, on introduit dans l'équation discrète (2 1 ) des solutions de la form e

r i, = ( r i • (4i)

9 S i le sch é m a à a n aly ser n 'est pas lin éaire à co efficients con stan ts, o n le rend te l e n a p p liquant les
s im p lific a tio n s q u e l'o n ju g e raiso n n ab le.

45
où v = / i A x . C om m e on considère q u e p < ( 2 k ) / ( 2 A x ) , on a
0 < V < n . (42)

L e facteur d'am plification des solutions d e ce type vaut

2 2 Af
y* = 1 — 2~ ^ “ cost>) • (4 3)
Ax

P o u r que le schém a soit stable, il faut |y*[ < 1, pour tou te valeur adm issible de v . D onc, les
incrém ents des variables indépendantes d oivent vérifier l'inégalité

^ 4 (44)
Ax 2 2
q ui est équivalente — en variables dim ensionnelles — à (38) quand K = Ax -1 -»«> .
O n va m aintenant utiliser largem ent l'analyse de von N eum ann p o u r appréhender les
p ropriétés essentielles d'autres schém as d e discrétisation de l'équation de la ch aleu r (15).
D ans (21b), la dérivée tem porelle de la tem pérature — le second m em bre de cette relation
— est évaluée à l'in stan t I A t, c'est-à-d ire l'in stan t le plus précoce co n cern é p a r l'éq u atio n
discrète. U n tel schém a est dit explicite. D ans un schém a im plicite, p ar con tre, on évalue la
dérivée tem porelle à l'in stan t le plus tardif concerné p a r l'équation discrète: p o u r transform er
(21) en u n schém a im plicite, il su ffit d'évaluer la dérivée tem porelle en (/+1) A t. Cette
discrétisation s'écrit
rp* rn* /T—
1* rji ♦ rp*

l+\,k ~ l l,k _ £ J/+l,k+l ~ l+ lk + l+\,k-\ (-4 5 ^)

At A x2
L e facteur d'am plification de ce schém a,

^ = - w i ----------: • (46)
1 + ------ (1 - cosd)
Ax
est toujours positif et inférieur ou égal à l'u n ité, quelles que soient les valeurs de la diffusivité
therm ique, de l'incrém en t tem porel et du pas spatial. En d 'autres term es, ce schém a e st toujours
stable. O n d it que ce schém a est inconditionnellem ent stable — alors que le p rem ier schém a était
conditionnellem ent stable.
L e p rix à payer p o u r bénéficier d 'u n e propriété aussi avantageuse que la stabilité
inconditionnelle est la nécessité de résoudre un systèm e d ’équations algébriques linéaires p o u r
p ro g resser d'un pas d e tem ps. E n ré-écriv an t l'équation (45) sous la form e

_ ^L p * + + - — - T* = T* (47)
Ax 2 { A x2 ) /+1’* A x 2 /+1’*+1 l>k ’ ( j

on v o it im m édiatem ent que la m atrice d u systèm e est tridiagonale. C 'est pourquoi il n ’e st pas
nécessaire d 'em p lo y er l'alg o rith m e d'élim ination de G au ss pour résoudre ce systèm e. On peut
se contenter de l'alg o rith m e de T h o m as, qui exige d'effectuer un nom bre d 'o p ératio n s
arithm étiques proportionnel à K , le nom bre de points d e discrétisation, alors que l'élim ination
de G au ss dem ande un n o m b re d'o p ératio n s proportionnel à i f 3. Le nom bre d 'o p ératio n s
exigées par le schém a im pücite étant du m êm e ordre de grandeur que p o u r le schém a explicite,
on p o u rrait conclu re qu'il est toujours souhaitable d'em ployer le schém a im plicite si sa précision
n 'était p as assez m odeste. En effet, cet algorithm e est d u prem ier ordre de précision selon le
tem ps et du second ordre selon l'esp ace — c'est-à-dire q u 'il n 'e st pas m eilleur, de ce point de
vue, que le schém a explicite.
Le schém a leapfrog — saute-m outon — ,
rrt*f r ji ♦ 'T ’ * r\ | r r '*

/+!.* ~ l-lk _ ^ Âl M l ~ l,k l l,k - 1 , .g,

2 At A x2 { '

estim ant la dérivée tem porelle de la tem pérature p ar une différence centrée, est d u second ordre
de précision dans le tem ps et l'espace. O n p eu t donc penser qu'il est supérieur aux deux
schém as en visagés p lus haut. M alheureusem ent, il est inconditionnellem ent instable, et donc
inutilisable! Voici un excellent exem ple du peu de pertinence des inform ations fou rn ies par
l'e rre u r de troncatu re...
L a solution la plus pratique pour augm enter la précision du schém a réside dans la
construction de schém as sem i-im plicites. En utilisant l'interpolation linéaire

Tl x . k “ ° “ X ) T '.t + X T M , k ' a v e c 0 S * £ >' <4S))

le schém a sem i-im plicite résultant de la com binaison d u schém a explicite (2 1 ) et d u schém a
im plicite (45) s'écrit
* * j* _ 2 j* , j*
U + i,k l i,k _ , l+ x M l i+ X'k - 1
 T ~ ~ A Ä ? ’ (50)
où x représente le niveau d'“im plicité” du schém a. Bien entendu, si x vaut zéro ou l'u n ité, on
retrouve le schém a explicite (21) ou le schém a implicite (45). L 'erreu r de troncature de
l'alg o rith m e (49) est du second ordre dans l'espace et d u prem ier ordre dans le tem p s, sauf
q u an d x = 1/2- D ans ce cas, l'erreur de troncature est du second ordre dans le tem ps. L e facteur
d'am p lification est
2 AAr
1 " (1 - 2 : ) 2" (1 - c o su )

= — — ü f - ---------- — - <51>
1 + x —T T C1 “ c o s'ü)
A x1
conduisant à la condition de stabilité
XAt < 1
Ax 2 ~ 2 m a x ( 0 ,l- 2 ^ )

A in si, dès q u e ^ > 0 , le schém a est “p lus stable” que le schém a explicite — c'est-à-d ire que le
p lu s grand p as de tem ps adm issible est p lus grand que celui associé au schém a explicite. Le
schém a est conditionnellem ent stable p o u r x < 1 / 2 , et inconditionnellem ent stable sinon —
c'est-à-dire quand la partie im plicite “dom ine” la partie explicite.
A u d elà de la stabilité, l'analyse de von N eum ann perm et de com parer le com portem ent du
facteu r d'am plification du schém a num érique, y * , avec celui de la solution exacte, y. O n obtient
ce dernier en introduisant dans l'équation de la chaleur continue (15) une solution du type
“m o d e de F o u rier”, c'est-à-dire

47
 T ( t ,x ) = e - £ ,* lfix , (53)

ce qui fournit
AA t
T ( t + A t,x ) = e - ta*-ï v
y = ---------------- , (.5. .4.)
T ( f,x )

^ A f/A x 2 = 0.01 A, A i/ Ax2 = 0.1

0 1/2 1

A .A Í/Ax2 = 1 A. Af / Ax2 = 10

0 1/2 1
D / 7C

Figure 7. C om paraison du facteur d'am plification de la solution continue (trait

plein épais), y, avec celui d e la solution discrète, y * , en fonction du no m b re d 'o n d e
norm alisé, i). On considère les schém as explicite (trait interrom pu), implicite
(pointillés) e t sem i-im plicite (trait plein) pour %= 1/ 2 .

D ans la F igure 7 , on v oit que le coefficient d'am plification de la solution discrète est
généralem ent très éloigné de celui d e la solution continue p o u r les m odes d o n t la longueur
d 'o n d e est courte, c'est-à-d ire proche de 2 A x . Cette constatation est valable p o u r quasim ent
to u s les schém as n u m ériq u es, q u 'ils v isen t ou non à résoudre l'équation de la ch aleu r. En
d 'au tres term es, les phénom ènes d ont la longueur caractéristique est proche de la taille de la
m aille so n t quasim ent tou jo u rs m al représentés par un schém a num érique. P ar c o n tre, tous les
algorithm es exam inés dans la F igure 7 — explicite, implicite et sem i-implicite — représentent
très bien l'évolution des m odes de grande longueur d 'o n d e. Ceci n 'est guère éto n n an t, ca r le
tem ps caractéristique d'am o rtissem en t de ces m odes est nettem ent p lu s grand que l'incrém ent
tem porel, A t, d e sorte que les détails de la discrétisation tem porelle im portent assez peu p o u r ces
m o d es.
U n e autre approche de la résolution num érique de l'équation de la chaleur est celle des
m éthodes d e m arches aléatoires, selon laquelle on considère un grand n om bre de particules
ponctuelles transpo rtan t chacune une quantité d'énergie interne appropriée avec une vitesse telle
que les d éplacem en ts de ces particules sim ulent la diffusion m oléculaire de la chaleur. Si X (t) est
la co o rd o n n ée spatiale d'une de ces particules à l'in stan t t, alors à l'in stan t t+ A t se trouvera à la
p osition

X (t + A t) = X ( t) + ( 6 AAt ) 1/2r , (5 5 )

où r est u n nom b re aléatoire, uniform ém ent distribué en tte - 1 et +1. L 'an aly se de cette
techn ique sort du cad re de la présente note.

E q u a tio n de la c h a le u r m u lti-d im e n sio n n e lle

D e no m b reu x p ro cessu s diffusifs ne peuvent être valablem ent m odélisés p a r une équation de
diffu sio n à une dim ension spatiale. On ne peut donc se limiter à l'étu d e d 'é q u atio n s et
d 'alg o rith m es uni-dim ensionnels: il faut considérer des problèm es où l'in co n n u e d ép en d d 'a u
m o in s deux coordo n n ées spatiales. A insi, il est souhaitable d'ex am in er une équation du type
dT . d2T , d2T
— = X — t + X — T , (56)
dt x dx2 y dy2 K }

où Xx et Xy représentent les diffusivités associées aux directions des axes des coordonnées
sp atiales x e t y. Si Ax e t A y sym bolisent les pas spatiaux selon les axes x et y , alors le schém a
num érique explicite le plus sim ple visant à la résolution de (56) est sans doute
rj-yfy m* A 'T’* J rTT->+
* l + l,kx ,ky ~ l,k :, k y U x+ \,k y ~ Z 1 l,kx ky + l l,kx - \ , k
_ ^ ------------------------------
At x A x¿

T* -2 T * +T*
l l,kx,ky + 1 V , i
+ ^ V ---------------- ’ (57)
où les indices en tiers k e t k servent à repérer la position des points de discrétisation selon les
x y
directions d es axes x et y . P ar une m éthode analogue à celle utilisée p o u r les p roblèm es uni-
dim ensionnels, on p e u t m ontrer que le facteur d'am plification de l'algorithm e (57) vaut
2X At 2 A Aí
y* = 1 —S
73— C1 -- co
“ (1 st) )
C0SV T-
7T ~ ^ _ c o s v y ) » (58)
Ax 2 * Ay

de sorte q u e la condition d e stabilité de ce schém a s'écrit

X At X At i

+ V s I ' (59)
S i les diffusivités et les pas spatiaux sont identiques dans les deux d irections, le p lu s grand

49
incrém ent tem porel adm issible est égal à la m oitié de celui découlant du critère (44), q u i prévaut
p our le schém a équivalent à une dim ension. Cette relative perte de stabilité s'accen tu e encore si
l'on co n sid ère trois dim ensions spatiales. C 'est pourquoi il est tentant d 'im p lém en ter un schém a
sem i-im plicite ou im plicite. C ependant, à deux dim ensions, un tel algorithm e exige la résolution
d 'un systèm e d'équations algébriques linéaires d ont la m atrice est pentadiagonale, et d ont deux
bandes so n t largem ent distantes de la diagonale principale. Com m e une telle opération peut
dem an der une q uantité prohibitive d'opérations arithm étiques, il p eu t être approprié d 'en v isag er
un algorithm e d'un autre type, à savoir la technique dite im plicite aux directions alternées. Cette
dernière consiste à p rogresser selon deux “dem i-pas de tem ps” avec un algorithm e implicite
dans une direction et explicite d an s l’au tre, la direction implicite changeant d 'u n dem i-pas de
tem p s à l'autre:
rp * rp * rp * r\ rp * rp *

1 l+ ll2 ,k x ,ky “ kx ,ky 1l+ l/2 >k x + l k y ~ l+ l/2,kxk y l M I 2 ,k x - U y

A t/2 * Ax 2

T l,k .k +1 ” 2 T l,k k + Tl , k . k - 1
+ X — ^ ( 6 0 a )
y Ay 2

T* -T * T* - 2 T* +T*
l+ \,k x ,ky 1 M I 2 .k x .ky _ ^ 1M /2 ,k x + hky l+ l/2,kx ky l M f2 ,k x - \,k y

A t/2 x A x2
* t_ ? T* + T*
l+ l,k ,k +1 l+ l,k k l+ i.k ,k -1
+ X ------- ^ ^ 4_>_ . (60b)
y A y2

Si le p as de tem ps est relativem ent grand, un des dem i-pas de tem ps peut être instable. M ais, en
to u t cas, le schém a com plet reste inconditionnellem ent stable. En effet, son facteur
d'am plification est

X At , \ Ar ,, ,
1 -- cost) ) 1 - - f y - (1 - co st) )
(1
y* _ Ax 2 ____________ A r __________ L ( 6 i)
7 “ X At X At ’ ^ ;
1 + - f r (! - co st) ) 1 + -2— (1 _ cost) )
Ax A y¿ y

une q uan tité d o n t la valeur absolue est toujours inférieure ou égale à l'unité.

E q u a tio n d 'a d v e c t io n -d iffu s io n

C om m e l'illustre la Figure 6 , la précision des schém as num ériques visant à fou rn ir une solution
approchée du problèm e de conduction considéré p lus haut est généralem ent b o n n e, sau f si le
schém a est instable. G rossièrem ent parlant, ici, il suffit quasim ent d 'assu re r la stabilité du
schém a p o u r obtenir une approxim ation num érique relativem ent p récise. Ceci est dû
principalem ent au fait que l'am o rtissem en t des m odes de F ourier de la solution est très fort.
D onc, il im porte finalem ent assez peu que l'am ortissem ent des m odes discrets soit ou ne soit
pas correct — p our autant q u 'ils restent stables. En fait, ce problèm e différentiel aux dérivées

50
partielles est san s do u te un des p lus simples à résoudre n um ériquem ent. Cette sim plicité cesse
d 'ê tre de m ise dès que l'o n considère des p rocessus p h y siq u e s qui ne cond u isen t p as à des
m o d es am ortis. A insi, si l'on passe d'un m ilieu solide à u n milieu fluide, il convient d 'a jo u te r à
l'équation de la chaleur (15) le transport convectif — ou “ a d v e c tif’. On o b tien t alors l'équation
dT_ _ _d_ dT
—u T + X —— I, (62)
dt d x ^ dx

où u sym bolise la vitesse du fluide. Si cette dernière est une constante p o sitiv e, on peut
adim ensionnaliser (62) de la façon suivante:
dT dT \ d2 T
(63)
dt dx Pe d x 2

où P e est un coefficient adim ensionnel généralem ent appelé “ nom bre de P éclet” . C e dernier
m esure l'im p o rtan ce du transport convectif p a r rap p o rt au transport d iffu sif — et a une
signification équivalente à celle du nom bre de R eynolds d a n s les équations de bilan de quantité
d e m o u v em en t d es fluides.

instant: t = 0 instant: t= 1
2

1
1 0 1 2 3

instant: r= 2 instant: t = 3
2 2
CD
•CD
Ui
ca
1 E 1 Í+TT
o
c
CD CD

3 0
•C
2?
•CD •CD
Q. CX
eCD E
CD

1 1
0 2 3 4 2 3 4 5
co ordonnée spatiale: x c o ord o n né e spatiale: x

F ig u r e 8. S olution exacte (trait plein) et solutions discrètes fou rn ies par

l'alg o rith m e ( 6 8 ) (o) et (69) (+) du problèm e d'advection uni-dim ensionnelle (64)-
(67) à différen ts instants. L 'instabilité du schém a centré ( 6 8 ) apparaît clairem ent,
to u t co m m e la dissipation erronée de (69).

51
 Si le transport a d ve c tif dom ine largement le transport d iffu s if, ob ten ir une approxim ation
num érique précise de l'équation (63) est très difficile, su rto u t si l'o n traite une v ersion multi-
dim ensionnelle de cette équation — c'est-à-dire si plus d u n e variable indépendante spatiale doit
être retenue. D es m illiers d'articles ont été consacrés à ce problèm e! O n est alors confronté à un
problèm e où il ne suffit plus — loin s'en faut! — d 'a ssu re r la stabilité de l'alg o rith m e. O n va
l'illustrer en exam inan t brièvem ent deux cas lim ites de l'équation d 'advection-diffusion (63).
Q uand le nom bre de Péclet est très g rand, il est tentant de négliger le term e diffusif. D ans
un dom aine infini, effectuer cette sim plification ne p ro v o q u e généralem ent p as de grandes
erreurs. C e faisant, on obtient l'équation hyperbolique

dt dx

d o n t on cherche la solution dans le dom aine -°« < x < +°°. C ette dernière est bien connue:
T ( t,x ) = T ( t = 0 , x - t ) . (65)

Le profil initial de tem pérature est donc transporté, san s déform ation, avec une vitesse
adim ensionnelle unitaire. Ici, on suppose que la tem pérature initiale est le “trait M orse” défini
par
T ( t = 0 ,x ) = 1, p o u r \x\ < 1 , ( 6 6 a)

7"(r = 0 ,x ) = 0 , p o u r 1 < |jk| . ( 6 6 b)

L a solution d e ce problèm e est alors

T ( t , x ) = 1, p ou r |x - t\ < 1 , (67a)

T ( t = 0 ,x ) = 0 , pour l < | x - r | . (67b)

On va m aintenant bâtir- deux approxim ations discrètes de la solution de (64), valables aux
in stan ts t[ = I A t (/ = 0, 1, 2, 3, ...) et aux points x k = { k - \ / 2 ) A x (k = 0, ± 1 , ± 2 , ± 3 ,...). En
utilisan t des techniques sem blables à celles em ployées p o u r le problèm e de diffusion de chaleur,
on p eut établir les schém as
m* rj~'* 'T* T*
I l+ l,k ~ 1 l ,k _ I l ,k + l ~ u - l /(TON

At 2Ax
rj-’^ rji♦ rr' * rriH1

l+ l,k ~ l,k 1 l,k ~ U - l

------ , (69)
At Ax
q ui d iffèrent par la façon d ont est discrétisée la dérivée spatiale. Le prem ier schém a contient une
approxim ation centrée de la dérivée spatiale, qui est précise au second ordre en À*.
M alheureusem ent, ce schém a est inconditionnellem ent instable, c'est-à-d ire q u 'il est instable
quelles q u e soient les valeurs de l'incrém ent tem porel et du pas spatial (Figure 8 ). Le second
schém a, basé sur une approxim ation décentrée — ou u p w in d — de la dérivée spatiale n 'e s t
précis qu'au p rem ier ordre en Ax, m ais est stable à condition que le nom bre de C ourant, A t / A x ,
so it inférieur à l'unité. L 'inconvénient m ajeur de ce schém a réside d an s le fort taux de

52
dissipation num érique qu'il contient, ce qui constitue u n e erreur significative 10 ca r la solution
exacte n 'est nullem en t dissipée ou amortie (Figure 8 ). P o u r obtenir une approxim ation discrète
acceptable d u problèm e d'advection ci-dessus, il faut a v o ir recours à des m éthodes relativem ent
so p h istiq u ées — d o n t la discussion sort du cadre de ces n o tes.
O n aborde m ain ten an t le second cas lim ite de l'éq u atio n d'advection-diffusion (63), à savoir
la situation stationnaire — c'est-à-dire indépendante du tem p s. L'équation à résoudre est alors
dT I d2T
- — + — - 0 . (70)
dx Pe d x

L e d o m ain e d 'in térêt est - 1 < x < 1 et les conditions aux lim ites sont
T (x = 0) = 1 , (71a)

T (x = l) = 0 . (71b)

L a solution du problèm e (70)-(71) vaut

1 _
T (x ) = — zpT ~ • (72)
~ e '

S i le nom bre de Péclet Pe est petit, alors le second term e de (70) tend à d om iner l'éq u atio n ,
qui tend à se réd u ire à ô 2 7’/3 x 2 = 0 , d ont la solution est T (x ) = 1 - x . P ar contre, si P e est grand,
le p rem ier term e de (70) tend à dom iner, de telle sorte q u e l'équation est à peu p rès d T / d x = 0.
L a solution d e cette dernière est T (x) = 1. A proxim ité de la frontière x = l , cette solution n 'e st
évidem m ent p as com patible avec la condition (71b) im posée à cette extrém ité du dom aine
d 'intérêt. L a so lu tio n s'ajuste alors vers la v aleur p rescrite à la frontière d an s une région de fort
gradient, appelée couche lim iten , dont l'ép aisseu r est de l'o rd re de P e~ l . L a F igure 9 illustre
clairem ent cette discussion.
L e p ro b lèm e discret com posé de l'équation
T -T i T - 2T + T
_ jk ± 1 l k - 1 + k+l 1 _ q «on
2 Ax Pe A x2

où T ,= T ( x = k A x ) — avec k = 0, 1, 2,... K = A x~ l — e t des co n d itio n s aux lim ites

7-0 = 1 , (74a)

T k = 0, (74b)

est com patible avec le problèm e différentiel (70)-(71). Il se fait que ce problèm e discret adm et,
p o u r P e * 2, une solution analytique simple:

10 L e sc h é m a u p w in d d u p rem ier ordre fo u rn it la so lu tio n exacte si A t / A x = 1. C eci e s t cependant san s grand

in té rê t, c ar le n o m b re d e C o u ran t d o it ê tre in férieu r à l’unité dès q u e la diffu sio n e st p ris e en c o m p te . E n o u tre , à
p lu sie u rs d im e n s io n s sp atiale s — c ’e st-à-d ire p o u r les p ro b lèm es réellem ent intéressan ts o u u tile s — , il n'existe
a u cu n e c o m b in a iso n d es in crém en ts des v ariab les indépendantes qui s o it su sc e p tib le d e c o n d u ire à la so lu tio n
ex ac te — s a u f d a n s d es c as très p articuliers.
11 L a n o tio n de “c o u ch e lim ite ” fu t san s d o u te développée p o u r la prem ière fo is en m écanique des flu id es, m ais
p e u t m a in te n a n t ê tre co n sid érée c o m m e u n co n cep t m athém atique a ssocié aux p ro b lè m e s de perturbation
s in g u liè re — que l'o n ren co n tre dans de n o m b reu x m odèles m athém atiques to ta le m e n t indépendants de la
m écan iq u e d e s fluides.

53
 K -k
2 - Pe Ax
1 -
2 + Pe Ax
k = 0 , 1 , 2 ,... K . (75)
Tk =
'2 - P e l * K
1
2 + Pe Ax

C ette expression présente un com portem ent oscillatoire s i Pe A x > 2. C om m e la solution exacte
n 'o scille p as, un tel com portem ent est difficilem ent adm issible. Il est donc souhaitable
d 'im p o ser q u e l'incrém ent spatial Ax satisfasse la co ndition

Ax < — . (76)
Pe
C om m e la lim ite supérieure adm issible — selon (76) — de l'in crém en t spatial est du m êm e
ordre de g ran d eu r q u e l'épaisseur de la couche lim ite, la condition (76) signifie q u e le pas spatial
ne p e u t être significativem ent plus grand que l'épaisseur d e la couche lim ite. En d 'au tres term es,
la couche lim ite d o it être résolue p ar la discrétisation choisie.

P e = 0.1 , P e A x = 0.01 P e = 1 , P e A x = 0.1

0.5 0 .5

AX = 0.1

P e = 1 0 , Pe Ax = 1 P e = 100 , P e A x = 10

0.5 0.5

c o o rd o n n é e sp a tia le: x c o o rd o n n é e sp a tia le : x

F ig u r e 9 . S olution exacte (trait plein) et solution discrète approchée (o) du

p ro b lèm e (7 0 )-(7 1) p o u r différentes valeurs du n om bre de Péclet Pe.

54
C o n c lu s io n

L es ex em ples ci-dessus illustrent les idées suivantes: .

— L a co nception d 'u n e m éthode de résolution n u m ériq u e d 'éq u atio n s différentielles aux
dérivées partielles n e s'im provise pas. Il faut m aîtriser u n certain no m b re de concepts de base,
com m e la com patibilité, la stabilité, la convergence. Il fa u t ég ale m en t être en m esure de choisir à
bon escient la fam ille du schém a à mettre en oeuvre, à sa v o ir le caractère explicite, implicite,
sem i-im plicite, etc.
— D 'au tre p art, p o u r traiter des problèm es réellem ent in téressan ts, c 'est-à-d ire des problèm es
beaucoup p lu s com plexes que les exemples académ iques abordés ici, il faut s'être forgé une
expérience su ffisan te. En effet, p o u r ces problèm es, au cu n résultat théorique rigoureux n 'e st
généralem ent disp o n ib le quant aux propriétés des sch ém as de résolution concevables. Par
co nséquent, seule l'expérience — assistée d 'u n bagage th éorique approprié — perm et de faire
les “ bons ch o ix ” d e m odélisation e t de discrétisation.
— L a sim ulation n u m érique d'un systèm e naturel ou artificiel dem ande la m aîtrise sim ultanée de
la “phy siq u e” du pro b lèm e, des m éthodes num ériques appropriées et de l'inform atique
nécessaire à la m ise en oeuvre sur ordinateur du schém a n u m érique choisi. E n o utre, il faut être
capable de réalise r la sym biose entre ces trois dom aines,
I l e st donc illusoire de p e n ser que la m ise en oeuvre d 'u n e m éthode num érique est une activité
de nature “technique", que l'on pourra apprendre au p ie d levé, le jo u r o ù le besoin s ’en fera
sentir. C ’e s t p o u rq u o i il est indispensable que le program m e de la licence en p h ysiq u e de nos
u niversités com porte une initiation aux m éthodes de sim ulation num érique de la p h y siq u e , avec
une p rio rité sensible p o u r la résolution de problèm es différentiels a u x d érivées partielles. Il faut
enseigner les no tio n s théoriques de base — les exem ples ci-dessus m ontrent d ’ailleurs q u ’il
s ’agit d ’une charge d ’enseignem ent plutôt m odérée — et, su rto u t, exiger que les étudiants se
form ent par la pratique, c ’est-à-dire en résolvant num ériquem ent, et de façon com plète, des
problèm es aux dérivées partielles.

R em erciem en ts. M ichel C rucifix a accepté de relire et co m m en ter la v ersion prélim in aire d e ce
texte. L ’au teu r a eu des discu ssio n s fructueuses à propos d u sujet abordé ici avec Jean-M arc
G érard, Ja n G o v aerts e t B ernard Piraux. L ’auteur est C hercheur qualifié au F onds N ational de
la R echerche S cientifique de Belgique.

R é fé r e n c e s

B ender C . M . and S. A. O rszag, 1978, A d va n ced Mathematical M eth o d s f o r S cien tists and
E n g in e e rs, M cG raw -H ill.
F erzig er J. H ., 1998 (2nd ed.), N um erica l M ethods f o r E ngineering A p p lica tio n , W iley.
H ack J. J., 1992, Clim ate system simulation: basic numerical & com putational co n cep ts, in:
C lim ate S ystem M o d elin g , K. E. Trenberth (E d.), C am bridge U niversity P ress, pp. 283-
31 8 .

55
H irsch C ., 1988, N um erical Computation o f Internal and E xternal F lo w s - V o lu m e 1:
F undam entals o f N um erical D iscretization, W iley.
R ichtm ëyer R'. D. and K . W . M orton, 1957, D ifference M ethods f o r Initial-V alue P ro b lem s,
W iley.
R o b in so n A . R ., 1987, Predicting open ocean cu rre n ts, fronts an d ed d ies, in: Three-
D im en sio n a l M odels o f M arine and E stuarine D y n a m ic s ,!. C. J. N ih o u l a n d B. M. Jam ait
(E d s.), E lsev ier, pp. 89-111.

56