CALCUL SCIENTIFIQUE 2 : DIFFRENCES FINIES POUR

QUATIONS AUX DRIVES PARTIELLES

PIERRE PUISEUX

1. Prliminaires

Un thorme.
Thorme 1. Le rayon spectral d'un matrice est un minorant de l'ensemble des normes
1.1.

matricielles subordonnes :

(A) 6 kAk
1.2.

EDP : notations et vocabulaire.

Cette section introduit des dnitions fonda-

mentales sur les quations aux drives partielles (EDP).

quation aux drives partielles : quation comprenant au moins une drive partielle
d'une fonction de plusieurs variables.

Ordre d'une quation aux drives partielles : ordre de la drive de plus haut degr
2
u
de l'EDP :
xu2 = f est d'ordre 2.
t
EDP linaire : quation de la forme

F (u, ux , uy , uxx , uxy , uyy , . . . , x, y) = 0

F est linaire par rapport aux variables, sauf ventuellement les deux
(x, y). Chacun des termes comporte au plus la fonction elle-mme ou une

o la fonction
dernires

de ses drives ; les coecients des termes peuvent dpendre des variables indpendantes. Exemples :

uxx + 3uxy + uyy + ux =uexy = 0

2
linaire, coecients variables : sin(xy)uxx + 3x uxy + uyy + ux =u = 0
2
x
non linaire : uxx + 3uxy + (ux ) =ue =y = 0

 EDP linaire, coecients constants :

 EDP
 EDP

1.2.1.

Classication sommaire des EDP linaires.

considrons une EDP linaire deux

variables, :

auxx + 2buxy + cuyy + dux + euy + f u = g

on lui associe son polynme caractristique en

et

P (, ) = a2 + 2b + c 2 + d + e + f
Les proprits de

P (, ) dterminent la nature de l'EDP Selon le discriminant b ac, le

polynme et l'EDP sont dits :

hyperbolique si b2 ac > 0
parabolique si b2 ac = 0
elliptique si b2 ac < 0

Exercice 2.

En TP, tracer les polynmes dans les trois cas. Par exemple

a=b=1
a = 1, b = c = 12
a = b = 1, c = 2
1

Janvier 2010

Figure 1.1.

Rsolution d'EDP par dirences nies

z = P (x, y)

Pierre Puiseux

de type hyperbolique, parabolique et elliptique

L'quation parabolique modle est :

1.2.2.

ut u = f

ut

dsigne la drive partielle par rapport au temps (u est donc une fonction de

variable d'espace, et de

t,

x,

variable de temps). Cette quation modlise par exemple la

conduction de la chaleur en rgime instationnaire.

L'quation elliptique modle est.

1.2.3.

u = f

(1.1)
o

u =12 u + 22 u , i ,

dsignant la drive partielle par rapport la i-me variable.

Cette quation modlise par exemple le phnomne de conduction de la chaleur stationnaire (c..d. en rgime permanent). En lasticit, on rencontre galement l'quation du
bi-laplacien, c..d. :

2 u = f

(1.2)

L'quation (1.1) peut tre discrtise par dirences nies, volumes nis ou lments nis. Les mthodes des dirences nies sont limites des domaines gomtriques simples .
L'quation (1.2) est le plus souvent discrtise par lments nis, pour des raisons de prcision.
1.2.4.

quation hyperbolique modle.

interviennent principalement en mcanique des uides

(aronautique, coulements diphasiques, modlisation de rupture de barrage et d'avalanches). Elles sont souvent obtenues en ngligeant les phnomnes de diusion (parce
qu'ils sont de faibles importance) dans les quations de conservation de la mcanique.
L'exemple le plus classique d'quation hyperbolique linaire est l'quation de transport
(ou d'advection).

(1.3)

ut ux = 0, t R, x R

avec condition initiale :

(1.4)

u (x, 0) = u0 (x)

Universit de Pau et des Pays de l'Adour

2 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Dans le cas o la condition initiale

u0

Pierre Puiseux

est susamment rgulire, il est facile de voir que

la fonction :

u (x, t) = u0 (x + t)

(1.5)

u0

est solution de (1.3)-(1.4). Si

est non rgulire (par exemple discontinue, il est possible

de montrer que la fonction dnie par (1.5) est solution en un sens que nous qualierons
de  faible .
1.2.5.

Conditions aux limites.

dlimite un domaine d'tude

Une EDP ne sut pas dnir une solution unique ; on

Rd

et on adjoint des conditions aux limites (C.L.)

l'EDP.

donne sur la frontire

u
Conditions de Neumann : la drive normale est donne
= a sur
n

u. n , o n est la normale extrieure .

Conditions aux limites Dirichlet :

Maillage.

u = u0

Avec

u
n

Rd . Un maillage de est un dcoupage de

en polydres (triangles ou en rectangles pour d = 2, ttradres, hexadres en dimension
d = 3) . Les sommets des polydres sont les nuds du maillage.

1.3.

Exemple.

on se

Soit un domaine born

= ]0, 1[ ]0, 1[,

donne n N et h =

1
. L'ensemble des points
n+1

h = {(ih, jh) , 1 6 i, j 6 n}
dnit un maillage cartsien rgulier (plus exactement le maillage est l'ensemble des
carrs dlimits par les droites d'quations

x = ih

et

y = jh,

avec

1 6 i, j 6 n)

Figure 1.2. Deux maillages cartsiens : irrgulier et rgulier, et un mail-

lage triangulaire

Universit de Pau et des Pays de l'Adour

3 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

On peut galement se donner deux suites de

(xi )16i6nx

et

(yj )16j6ny .

[0, 1],

Pierre Puiseux

nies et strictement croissantes

On considre alors le maillage cartsien irrgulier

h = {(xi , yj ) 1 6 i 6 nx , 1 6 j 6 ny }
Plus exactement le maillage est l'ensemble des rectangles dlimits par les droites

x = xi

et

y = yj , avec 1 6 i 6 nx , 1 6 j 6 ny ). Ponctuellement, on peut avoir besoin

x0 = y0 = 0 et xnx +1 = yny +1 = 1
une EDP linaire avec ses CL sur un domaine

d'ajouter des points de la frontire en posant,

tant donne

Au = b
usuellement on ne sait pas la rsoudre analytiquement. Dans ce cours, on se contentera
de trouver, par la mthode des dirences nies, une approximation de la solution aux
N
points du maillage h . Cette approximation, note Uh est donc un vecteur de R
o N
est le nombre de points du maillage.
Si une des variables est le temps, on la note

t [0, T ]

(plutt que

ou

y ),

et de la

mme manire que pour les variables d'espace, on discrtise le temps en considrant une

tn [0, T ], que l'on prendra le plus souvent quirpartis, c'est

tn = n.t o t > 0 est le pas de temps.
conditions, si u : [0, T ] R est une fonction, la valeur u (tn , xi , yj ) est

suite croissante d'instants

dire de la forme
Dans ces
uni,j .

note
1.4.

Les principales formules de dirences nies.

On rappel le thorme de Taylor

sous la forme de Young :

Thorme 3. soit f une fonction de classe

C n+1 au voisinage d'un point x R Alors

pour tout h rel,

f (x + h) =

(1.6)

n
X
hk

k!

k=0

f (k) (x) + O hn+1

La formule (1.6) est appele formule de Taylor d'ordre k (pour f au voisinage de x)

(k)
Ici f
(x) dsigne la drive k i`eme de f, et O (hn+1 ) est une quantit qui tend vers 0
n+1
aussi vite que h
 lorsque h tend vers 0.


 (h) 
n+1
En toute rigueur, on dit que (h) est un O (h
) si et seulement si  hn+1  reste born
au voisinage de

h = 0.

La formule de Taylor donne une approximation de

drives au point

f (x + h) ,

connaissant

et ses

x.

Rciproquement, connaissant la fonction

approximation des drives successives de

au voisinage de

au point

x,

on peut en dduire une

x.

Proposition 4. Les formules suivantes donnent des approximations des drives de f au

point x :

(1.7)

(1.8)

(1.9)

(1.10)

f (x + h) f (x)
+ O (h)
h
f (x) f (x h)
+ O (h)
f 0 (x) =
h


f (x + h) f (x h)
f 0 (x) =
+ O h2
2h


f (x + h) 2f (x) + f (x h)
2
f 00 (x) =
+
O
h
h2
f 0 (x) =

Ces formules sont la base des mthodes de dirences nies utilises pour rsoudre les
quations aux drives partielles de la physique. Elles se dmontrent en crivant la formule
de Taylor un ordre ad hoc, aux points
Universit de Pau et des Pays de l'Adour

xh

et

x + h,

et en faisant des combinaisons

4 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Pierre Puiseux

linaires de ces quations. Pour la troisime formule, par exemple, on utilise

points

x+h

et

xh

(1.6)

aux

h2 00
f (x) + O h3
2


h2
f (x h) = f (x) hf 0 (x) + f 00 (x) + O h3
2
f (x + h) = f (x) + hf 0 (x) +

(1.11)

(1.12)

en soustrayant la seconde quation de la premire , on obtient

2hf 0 (x) = f (x + h) f (x h) + O h3
O(h3 )
h

d'o, en remarquant que

Exercice 5.
Si
avec

= O (h2 ) ,

on dduit la formule anonce.

Dmontrer les autres formules proposes.

f
2f
Rd R alors les formules s'appliquent aux drives partielles x
et
xi xj
i

f :
1 6 i, j 6 d

Exemple 6.

Pour

d=2

on a :

2f
1
2
(x,
y)
=
(f
(x
+
h,
y)

2f
(x,
y)
+
f
(x

h,
y))
+
O
h
x2
h2
2


1
f
(x, y) = 2 (f (x, y + h) 2f (x, y) + f (x, y h)) + O h2
2
y
h
on en dduit que

f (x, y) =

1
2
(f
(x
+
h,
y)

4f
(x,
y)
+
f
(x

h,
y)
+
f
(x,
y
+
h)
+
f
(x,
y

h))
+
O
h
h2

ce qui se traduit plus visuellement en disant que la molcule (stencil) du laplacien de

pas

est

1
1
= 2 1 4 1
h
1

h
1.5.

Exemple basique.

Dans le cas gnral, on cherche rsoudre une EDP avec ses

conditions aux limites

Au = b
l'oprateur

prend en compte les conditions aux limites. Par exemple

On considre l'EDP

(1.13)

00
0

u + cu
u (0)

u (1)

=f
=0
=0

dans

f, c

]0, 1[

sont des fonctions donnes,

1
de longueur h =
n+1

c > 0.

On dcoupe le segment

[0, 1]

en

n+1

intervalles

h = {ih, 0 6 i 6 n + 1}

Les conditions aux limites s'crivent :

u0 = un+1 = 0

L'approximation par dirences nies de la drive seconde au point

00

ui

ih

est

1
2
(u

2u
+
u
)
+
O
h
i1
i
i+1
h2

1. En ralit c'est une EDP une seule variable, donc une simple quation direntielle, mais elle
permet de bien comprendre le fonctionnement de la mthode des dirences nies.
Universit de Pau et des Pays de l'Adour

5 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Pour la drive premire, on prend une approximation en

Pierre Puiseux

O (h2 )

an de garder la

mme prcision que sur la drive seconde

u0i =

1
(ui+1 ui1 ) + O h2
2h

Avec ces trois quations, on peut rcrire (1.13) aux points du maillage :

1
c1
2
(u

(u

2u
+
u
)
+
u
)
=
f
+
O
h
2
0
1
2
0
0
h2
2h


1
ci
(1.14)
1 < i < n : 2 (ui1 2ui + ui+1 ) +
(ui+1 ui1 ) = fi + O h2
h
2h


1
cn
i = n : 2 (un1 2un + un+1 ) +
(un+1 un1 ) = fn + O h2
h
2h
Ce qui peut s'crire matriciellement Ah Uh = Bh , soit, en extension :

c1
2
h12 + 2h
0

0
u1
h2




c2
c2
2
1
h12 2h

0
+
2
2

h
h
2h

u2
.
..
..

..
.
.
.
0
.

.

..
cn1
1
2
.

.
h2 + 2h un1
.
h2


cn
2
un
0

h12 2h
h2
i=1 :

La rsolution de ce systme fournit

Uh = A1
h Bh

f1

f2

= ...

qui contient les approximations de la

u aux points du maillage.

k
k = =2 + 2 cos n+1
, k = 1, 2, . . . , n;

solution

1.6.

Consistance, stabilit, convergence.

La premire question qui se pose concerne

la qualit de la solution : est-elle proche de la solution exacte ?

0 (ou

bien n vers l'inni), la solution approche Uh converge-t-elle vers les valeurs exactes Uh =

(u (ih))16i6n de la solution ?
La seconde question concerne la convergence : si l'on fait tendre le pas

vers

De manire gnrale, on considre une EDP avec ses conditions aux limites et initiales :

(u) = 0

(1.15)
de solution

u : 7 R

et

h ,

un maillage de

gnralement dans ce cours, le domaine

est :

= ]0, 1[ ]0, +[

(x, t) , x est la variable d'espace et t la

variable temps. Le maillage est de la forme h = {(jx, nt) , 1 6 j 6 J , et n > 0}.
On a pos h = (x, t)
ou bien = ]0, 1[ ]0, 1[ si le temps n'intervient pas et le maillage est de la forme
h = {(ix, jy) , 1 6 j 6 J , et 1 6 i 6 I}. On a pos h = (x, y)
En chaque point (jx, nt) l'quation (u) = 0 est discrtise par dirences nies sous
, o pour tout point

une forme

nj (u) = 0

(1.16)

Dnition 7.

Troncature

On appelle erreur locale de consistance ou erreur locale de troncature au point

(jx, nt)

la valeur :

nj = nj (
u)
o

et la solution exacte

L'erreur (globale) de troncature est



(x, t) = max nj , (jx, nt) h
Universit de Pau et des Pays de l'Adour

6 Laboratoire de Mathmatiques Appliques

n1

fn

Janvier 2010

Rsolution d'EDP par dirences nies

Pierre Puiseux

Le schma est consistant si

lim

(x, t) = 0

(x,t)(0,0)

Le schma est consistant d'ordre

en espace,

en temps si

(x, t) = O (xp ) + O (tq )

Dnition 8. Si les fonctions nj sont linaires, alors la mthode peuts'crire


comme un
systme linaire
indpendament

Dnition 9.

Ah Uh = Bh .
de h

La mthode est stable si et seulement si

On appelle restriction la grille

A1 
h

est born

l'application

h : C 1 (]0, 1[) Rn
u 7 (u1 , u2 , . . . , un )t
et

h (u)

sera gnralement not

Uh

ou parfois mme

uh .

Thorme 10. Une mthode stable et consistante est convergente.

Dmonstration.

Ah Eh

Eh = A1
h Rh

donc

Eh = Uh Uh
= Ah Uh Ah Uh
= Ah Uh Bh
= Rh

TODO A revoir...L'erreur

vrie

puis



 kRh k
kEh k 6 A1
h

2. Diffrences finies pour les problmes elliptiques
2.1.

En dimension 1.

(2.1)

2.1.1.

00

u = f
u (0) = 0

u (1) = 0

sur

On reprend la problme modle du premier chapitre :

]0, 1[

Discrtisation. h =

prcdent, on pose

1
,n
n+1

> 1,

on dnit la discrtisation comme au chapitre

u (ih) = ui

La discrtisation de (2.1) amne le systme linaire

Ah Uh = Bh

(2.2)
avec

Ah =
Uh = (ui )16i6n et Bh = (fi )16i6n
e`me
La i
composante de l'erreur

1
tridiagn (1, 2, 1)
h2

1
(vi+1 + 2vi vi1 ) vi00 . Or,
h2
on sait que (en utilisant la formule de Taylor-Lagrange) pour tout i, 1 6 i 6 n on a
de troncature est

ri = h2

ri =

v (4) (i ) + v (4) (i )
24

Universit de Pau et des Pays de l'Adour

7 Laboratoire de Mathmatiques Appliques

Janvier 2010

avec

constante

Rsolution d'EDP par dirences nies

](i 1) h, (i + 1) h[. Par

indpendante de h, on peut crire
et

dans

suite, en posant

Pierre Puiseux

C =

1
12

 
sup]0,1[ v (4) 

Rh = h2 (K1 , K2 , . . . , Kn )t

(2.3)

|Ki | 6 C, 1 6 i 6 n

Analyse de convergence en norme euclidienne

k.k2 La norme matricielle subordonne la norme euclidienne est la norme spectrale, qui, dans le cas d'une matrice A
2.1.2.

symtrique est gal son rayon spectral

(A) = max {|| , sp (A)}

o
n


1
1
1
La norme euclidienne de Ah est donc Ah
= max || , sp (Ah ) =
c'est dire

sin
h
2 (n + 1)

2
indpendamment de h

A1
h

donc

1
min{||,sp(Ah )}

est borne

A1
h

L'erreur de troncature est

kAh Vh (Av)h k2 = h

sX

Ki2

16j6n

6 h nC

6 h hC
2

La mthode est donc convergente en norme euclidienne.

2.1.3.

Analyse de convergence en norme innie k.k .

Dans ce cas, l'erreur de troncature

est

La norme matricielle est

kAk

kAh Vh (Av)h k 6 h2 C
P
= maxi j |aij |. On montre au paragraphe suivant (2.2)que
1
A 6 1
h
8

La mthode est donc convergente en norme innie.

Monotonie, stabilit en norme k.k .

Dnition 11. Soit A = (ai,j ) Rn,n .

2.2.

On dit que
On dit que
On dit que

A
A
A

Proposition 12.

ai,j > 0, i, j {1, . . . , n}

1
est monotone si A est inversible et A
>0
n
conserve la positivit si v R , Av > 0 = v > 0
est positive, et on note

A>0

si

A est monotone A conserve la positivit

Dmonstration.

Si A conserve la positivit, alors. Supposons Ax = 0, alors Ax > 0 donc

x > 0. de mme Ax 6 0 donc x 6 0 donc x = 0 donc A est inversible. En posant y = Av ,
1
e`me
1
1
on a y > 0 = v = A y > 0. Or la i
colonne de A
est A ei > 0, o ei > 0 est le
e`me
1
i
vecteur de base. Donc A
est positive et A est monotone.
1
Rciproquement si A est monotone, alors les coecients de A
sont positifs, donc si
1
y = Ax > 0 alors x = A y > 0, A conserve la monotonie.


Proposition 13. La matrice Ah = h1 tridiag (1, 2, 1) est monotone.

2

Universit de Pau et des Pays de l'Adour

8 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Dmonstration.
indice tel que

v Rn tel que Ah v > 0.

vp = min16i6n vi . on a donc
Soit

On pose

v0 = vn+1 = 0.

Pierre Puiseux

Soit

le plus petit

vp1 > vp
vp+1 > vp

(2.4)
(2.5)
On suppose

vp < 0.

Ah v > 0,

Comme

on a

(vp vp1 ) + (vp vp+1 ) > 0

ce qui est en

contradiction avec (2.4) et (2.5)

Proposition 14.

 1 
A  6
h

Dmonstration. A1
> 0
h
t

eh = (1, . . . , 1)

1
8

car

donc la mthode est stable en norme k.k

Ah

bi,j = A1
h

est monotone. De plus, en posant

i,j

et

1
, la norme de Ah s'crit :

 1 
A 
h

= max

bi,j = A1
h eh

x(1x)
est solution de (2.1) avec e = 1 au second membre. Comme
Or on voit que u (x) =
2
(4)
u = 0, l'erreur de troncature est nulle, donc uh = h (u) = (u (h) , u (2h) , . . . , u (nh))t
1
vrie exactement Ah uh = eh , donc uh = Ah eh et on obtient l'estimation

 1 
A 
h

6 sup |u|
]0,1[

1
8


Inuence de la discrtisation des CL.

2.3.

(2.6)

00

u = f
u0 (0) =

u (1) = 0

2.3.1.

Consistance.

i
0
1
16i6n
.
.
.

sur

Considrons :

]0, 1[

On discrtise :

Ah uh
(Au)h
Rh
0
(u1 u0 )
u (0) =
hK0
1
u (h) = f1 h2 K1
(u0 + 2u1 u2 )
h2
1
(ui1 + 2ui ui+1 ) u (ih) = fi h2 Ki
h2
1
h

.
.
.

.
.
.

n+1
un+1
1
Avec K0 =
u (h0 ) , Ki =
2
pour 1 6 i 6 n

h2
24

.
.
.

un+1 = 0
0
 (4)

(4)
u ((i i ) h) + u ((i + i0 ) h)

et

0 , i , i0 [0, 1]
t

Rh = (hK0 , h2 K1 , . . . ,h2 Ki , . . . , h2 Kn , 0)
ou les

1
h par K = max 12
sup[0,1] u(4)  , 12 sup[0,1] |u| .

Ce qui donne une erreur de troncature

|Ki |

sont borns indpendament de

Pour rsoudre le systme, on teste plusieurs possibilits :

On conserve la variable

u0 et on limine la variable un .

Le systme s'crit alors :

h
Ah Uh = B

avec

h h
1 2 1

1
..
..

.
.
Ah = 2
h
..

.
avec

h = (, f1 , f2 , . . . , fn )t
B

et

.
Rn+1,n+1

..
. 1
1 2
..

Uh = (ui )06i6n .

Universit de Pau et des Pays de l'Adour

9 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Le problme est que probablement

1
Ah

Pierre Puiseux

n'est pas borne. Donc le shma n'est pas

convergent.

1
, sous la
On peut changer la premire quation de ce systme. en la multipliant par
h
1

forme 2 (u1 u0 ) =
. Dans ce cas, le systme s'crit Ah Uh = B h avec
h
h

Ah

1 1
1 2 1

1 2

Rn+1,n+1

t

, f1 , f2 , . . . , fn et Uh = (ui )06i6n . ce qui revient crire la premire quation

avec B h =
h
1
(u1 u0 ) = h + hh K0 . Dans ce cas, le conditionnement de la matrice est meilleur, mais
h2
l'erreur de troncature devient en K0 c'est dire que le schma n'est pas consistant.
(L'erreur de troncature ne tend pas vers 0). Le shma n'est donc toujours pas convergent.
Le conditionnement
des matrices

2 1
1 2 1

Ah =

1 2

1 1
1 2 1

Ah =

et

Ah

est le suivant :

1 2

h h
1 2 1

Ah =

Ah , Ah

1 2

Universit de Pau et des Pays de l'Adour

10 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Pierre Puiseux

4133.6429 16370.242 235490.39

16373.242 65166.524 1319338.8
36718.536 146385.58 3623825.2
65169.524 260027.42 7426799.1
101726.21 406092.04 12961285.
146388.58 584579.44 20432189.
199156.66 795489.62 30024595.
260030.42 1038822.6 41908618.
329009.89 1314578.3 56242555.
406095.04 1622756.8 73175050.
On peut montrer, malgr cette inconsistance (ou cette instabilit suivant le cas), que
globalement, le shma reste quand mme convergent d'ordre 1, en norme

Une autre possibilit est de discrtiser la condition aux limites de Neuman l'ordre

h2
u (0) + O h3
2
2


h
= h + f0 + O h3
2

u (h) u (0) = hu0 (0) +

u1 u0
1
h2

+ 12 f0 + O (h). Le systme linaire s'crit dans ce cas :

h

t
f0
,
f
,
f
,
.
.
.
,
f
. On a alors la consistance en O (h) et la
1
2
n
2
matrice moins mal conditionne.

(u1 u0 ) =
h avec B
h = +
Ah Uh = B
h
et nalement

TP :.

Rsoudre le problme l'ordre 1 puis 2

00
2

u = sin x
u0 (0) =

u (1) = 0

(2.7)

sur

]0, 1[

u = sin x
1
1
, assemble la matrice Ah = 2 Tridiagn (1, 2, 1)
Ecrire une fonction qui, pour n donn, h =
n+1
h
et la matrice Ah du pb ci dessus. Calculer les spectres de Ah et
de Ah et tracer sur
1
1




un mme graphique les normes kAh k2 , Ah
, Ah 2 et Ah
2
Comparer avec la solution exacte

2.3.2.

Stabilit, tude de Ah .

On ne connat pas les valeurs propres de

Ah ,

on sait par

Rh = Bh Ah Uh l'quation

contre, tudier sa norme

Posons

2.4.

En dimension 2.

Assembler matrice et trouver valeurs propres.

Numrotation des inconnues

Thorme 15. Les valeurs propres de

Ah laplacien 2d de stencil

sont
pq
(p,q)

ui,j

1
h2

1
1 4 1
1



2 p
2 q
sin
+ sin
n+1
n+1
ip
jq
= sin
sin
n+1
n+1

4
= 2
h

Mettre en vidence la structure par blocs

Universit de Pau et des Pays de l'Adour

11 Laboratoire de Mathmatiques Appliques

Janvier 2010

2.5.

Rsolution d'EDP par dirences nies

Pierre Puiseux

Exercices.

Exercice 16.
monotone. Soit

Il s'agit de montrer que la matrice

v Rn

1 1
1 2 1

Ah =

1 2

est

tel que

Ah v > 0

(2.8)
Considrons le plus petit indice

p, 1 6 p 6 n

tel que

vp = min16i6n vi

et supposons que

vp < 0

(2.9)
(1) Montrer que

p=1

(2) On suppose que

conduit une absurdit.

p > 1.

tablir que

vp1 > vp
vp+1 > vp

(2.10)
(2.11)

et dduire de (2.8) une absurdit.

(3) En dduire que

Exercice.

Ah

A = tridiagn (1, 2, 1).

ik
ui = sin n+1
, 1 6 i 6 n.

Etude de la matrice

le vecteur de coordonnes

est monotone.

Montrer que

Soit

est un vecteur propre de la matrice

A.

x

16k6n

et

u Rn

Calculer la valeur propre

associe.

En dduire que

est inversible, calculer la norme spectrale de

conditionnement de

celle de

A1 ,

et le

en norme spectrale.

Etudier le comportement de ces direntes quantits lorsque

Exercice.

A,

n .

kAk :
Au 6 0 = u 6 0 (cette ingalit s'entend composante par composante. Considrer un indice i tel que ui = max16j6n uj et supposer ui > 0. En
distinguant les cas i = 1,1 < i < n et i = n parvenir une contradiction.)
En dduire que A est inversible (Au = 0 = u = 0)
Montrer que (A1 )ij = bij > 0 (considrer ei = A (A1 ei ) = Aci > 0 o ei est le ie`me
e`me
1
vecteur de base et ci est la i
colonne de A .
t
1
Montrer que kA1 k 6 8 . (Considrer u (x) = 21 x (1 x) etu = (u (h) , u (2h) , . . . , u (nh))
,

P
P
P
t
1
montrer que Au = (1, 1, . . . , 1) U. Comme u = A U =
j b1j ,
j b2j , . . . ,
j bnj ,

Etude de

Montrer que

on a

1
A
= kuk

6 inf u
[0,1]

1
8

3. Diffrences finies pour les problmes hyperboliques

pas de conditions aux limites (

= R),

Universit de Pau et des Pays de l'Adour

condition initiale
12 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Pierre Puiseux

Equation de convection, pb modle.

3.1.

u
u
+c
= 0 pour (x, t) R R+
t
x
u (x, 0) = u0 (x) donn

(3.1)
(3.2)

Thorme 17. Si

u0 est de classe C 1 sur R, alors l'quation

unique u (x, t) = u (x ct)

(3.1)

admet la solution

Dmonstration.

La fonction

0 0

0 0

u (x ct) vrie bien ut + cux = c (u ) + c (u ) = 0 et u (x, 0) = u0 (x).

L'unicit se dmontre en utilisant les droites caractristiques : on appelle droite

x0 ,la

caractristique de (3.1), issue de

droite

Dx0

du plan

(x, t),

d'quation

x (t) =

x0 + ct. Montrons que si u est solution, de (3.1), alors u est constante le long de cette
(x, t) Dx0 ), et ceci pour tout x0 . La fonction : t 7 u (x0 + ct, t)

droite (i.e. pour

vrie :

0 = cux + ut
= 0
donc

est constante et vaut

(t) = (0)

c'est dire :

(x0 , t) R R+ , u (x0 + ct, t) = u0 (x0 )

et en posant
3.2.

x = x0 + ct,

on obtient

u (x, t) = u0 (x ct).

Stabilit L2 , analyse de Von-Neumann.

En gnral, pour la stabilit d'un schma,

on utilisera l'analyse de Von Neumann que l'on dcrit ci dessous :

soit

est une fonction

2 -priodique

susament rgulire, on peut l'crire comme une

ikx

srie de Fourier, somme dnombrable d'ondes monochromatiques

u (x) =

u (k) eikx

kZ
avec des coecients de Fourier en nombre inni, dnombrable

u (k) =
(
u (k))kZ

1
2

kZ

u (x) eikx dx

est une suite double de coecients de Fourier.

De manire analogue, si

est intgrable sur

R,

elle est la somme d'un nombre inni

non dnombrable d'ondes monochromatiques

u () eix d

u (x) =
R

avec les coecients de Fourier cette fois ci en nombre inni, non dnombrable (c'est
dire une fonction

appele la transforme de Fourier) :

u () =

1
2

u (x) eix dx

R
L'application

u 7 u

est appele transformation de Fourier, la fonction

transforme de Fourier de
On retiendra que

est appele la

u.

u est un signal compos de frquences (ou k ), et que u () (ou u (k))

(ou la k -ime frquence) dans

reprsentent l'intensit, ou l'importance de la frquence

la fonction

An d'tudier le comportement d'un schma comme (3.6), on l'crit sous la forme
(3.3)

u(n+1) (x) = u(n) (x)

(n)
u (x + x) u(n) (x x)
2

Universit de Pau et des Pays de l'Adour

13 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Pierre Puiseux

x R= jx. Supposant que l'on peut dcomposer un , un+1 ,etc. en intgrale de Fourier,
un (x) = R un () eix d , on peut se demander, pour une frquence donne, que devient
le coecient de Fourier un () au cours de la marche en temps du schma ?
On dmontrera en exercice que la fonction : x 7 u (x + x) a pour transforme de
b () = u () eix
Fourier
avec

Lorsque l'on prend la transforme de Fourier de (3.3) (c'est dire lorsqu'on calcule le

), on obtient :




ix
ix
d
\
(n+1)
(n)
u
() = u () 1
e
e
2
(n) () (1 i sin x)
= ud

coecient de Fourier correspondant la frquence

La quantit

A () = 1 i sin ix

est appelle coecient d'amplication du schma.

On a donc





\
(n+1) () = u
0 () An+1 ()
b
u





n+1



 b0
6 u () sup |A|
R

Lorsque ce coecient est de module plus petit que 1, pour toute valeur de
L2 stable.

le shma

est

Lorsqu'il existe des frquences

pour lesquelles

A () > 1,

le schma est

L2 instable.

C'est le cas dans l'exemple ci-dessus.

2
On parle de stabilit L car la transforme de Fourier est une isomtrie et vrie :

2
(n) 2

d
(n)
u 2 =
u 2
L
Z  L
2

d
(n)
=
u () d
ZR

2

2n  d
(0)
=
|A ()| u () d
R

2
d
(0)
6 sup |A ()|2n u
2
L

On a donc

(n)
u 2 6 kAkn u(0) 2

L
L
3.3.

Amplication, dissipation, dispersion TODO.

L'analyse de Von Neumann donne

plus de renseignements que l'amplication globale du schma. Elle renseigne galement

sur la dissipation et la dispersion. La dissipation caractrise la facult d'un schma dissiper l'nergie d'une onde de frquence donne, la dispersion rend compte de l'aptitude de
ce mme schma disperser au quatre vents les direntes harmoniques de l'onde initiale
qui ne se propagent pas toutes la mme vitesse.
La solution exacte de l'quation de transport au temps

nt,

s'crit

un (x) = u0 (x cnt)
sa transforme de Fourier (en espace) de frquence

est

un () = u0 () eict

n

Le nombre

() = x
Universit de Pau et des Pays de l'Adour

14 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

est appel la phase de la frquence

Pierre Puiseux

dans le schma. C'est un nombre que connaissent

t
bien les physiciens. On le fait apparatre dans l'quation compte tenu de = c
:
x

un () = u0 () ei()

(3.4)

n

Lorsque l'on calcule la transforme de Fourier de la solution numrique au temps

nt,

on trouve

un () = u0 () An ()

(3.5)
En mettant

A ()

sous la forme

A () = |A ()| ei()()
i()
on en vient naturellement comparer les deux coecients e
et

quotient
est l'erreur de dispersion ... TODO, cf quarteroni, page480

3.4.

Shma explicite centr.

(n+1)

uj

(3.6)

on pose

Consistance en

ei()()

: le

t
= c x

(n)

= uj


 (n)
(n)
uj+1 uj1
2

O (t) + O (x2 )

2
Le schma est inconditionnellement L instable car : on crit (3.6) sous la forme

un+1 (x) = un (x) 2 (un (x + x) un (x x)) et on calcule sa transforme de

Fourier on obtient :





ix
ix
n+1
n
d
c
e
e
u
() = u () 1
2
d'o l'on tire :

n+1 () = A () u
cn ()
ud
avec le coecient d'amplication

1. Le schma est
2 sin2 x > 1
3.5.

inconditionellement stable.|A ()|

3.6.

1
1+ 2 sin2 x
2

= |1 i sin ix| = 1 +

Shma implicite dcentr amont.

n+1
un+1
= unj un+1
j
j+1 uj

(3.7)

A () = 1 i sin x et |A ()|2 =

Erreur de troncature : localement d'ordre 1 en temps, 1 en espace.

2
Le schma est inconditionnellement stable L car le coecient d'amplication vaut

2
 1 
1
= 1 + 2 (1 cos x) + 2 2 (1 cos x) > 1
A () = 1+(eix
et 
1)
A() 

Stabilit l
Principe du maximum
Variation totale

Shma avec diusion numrique articielle.

(3.8)

3.6.1.

(n+1)

uj

(n)

= uj

Erreur de troncature




 (n)
(n)
(n)
(n)
(n)
uj+1 uj1 + d uj+1 2uj + uj1
2

Voir le schma de Lax-Wendro

Universit de Pau et des Pays de l'Adour

15 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Stabilit L2 :

Pierre Puiseux

t
A () = 1ic x
sin x2d (1 cos x).
Pour c et d xs, A () dcrit dans le plan complexe une ellipse paramtre par [0, 2].
Les axes de cette ellipse sont vertical et horizontal. Le max de |A ()| est atteint pour


. Le schma est stable si et seulement si
x 0, 2 , , 3
2

3.6.2.

Le coecient d'amplication est

|1 4d| 6 1
|1 i 2d| 6 1
si la condition CFL est ralise :

x
> c
t
on obtient la condition de stabilit

1 1p
1
>d>
1 2
2
2 2
3.6.3.

Stabilit l .

On crit le schma sous la forme






(n)
(n)
(n)
uj1 + (1 2d) uj + d
uj+1
=
d+
2
2

(n+1)
uj

(n)

(n)

(n)

= Auj1 + Buj + Cuj+1

Les conditions obtenues pour

B > 0.

ci-dessus impliquent que

Si on suppose de

surcrot

d >
alors

A, B

et

,
2

sont positifs si bien que









 (n) 
 (n) 
 (n) 
 (n+1) 
uj
 6 A uj1  + B uj  + C uj+1 


 (n) 
6 (A + B + C) max uj 
j

ce qui assure la stabilit au sens

3.7.

Schma saute-mouton.
(n+1)
uj

(3.9)

l .

(n1)
uj

(n)
uj+1

(n)
uj1

consistance : le schma est d'ordre 2 en temps et 2 en espace car on reconnat des

drives centres en temps et en espace, d'ordre 2 donc :
(3.9)

Stabilit

L2



1  (n+1)
c  (n+1)
(n1)
(n1)
u
uj
+
u
uj
=0
2t j
2x j
 (n)
 (n)




u
u
+c
+ O x2 + O t2 = 0
t j
x j

: on crit (3.9) sous la forme suivante :

u(n+1) (x) = u(n1) (x) u(n) (x + x) u(n) (x x)
avec

x = jx.

En prenant la transforme de Fourier :

d
\
(n+1) () = u
(n1) () 2i u
(n) () sin x
u\
Cette galit couvre deux pas de temps. On peut l'crire matriciellement :

(n+1) ()
u\
(n) ()
ud


=

2i sin x 1
1
0

Universit de Pau et des Pays de l'Adour

(n) ()
ud
(n1) ()
u\

16 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Pierre Puiseux

Le coecient d'amplication apparat comme une matrice d'amplication


A () =

2i sin x 1
1
0

det (A () I) = 2 + 2i sin x
n 1. Sous la p
> c, le discriminant est positif et les deux racines 1 , 2 i sin x 1 2 sin

dont les valeurs propres sont racines de

condition CFL
vrient
3.8.

x
t

|1 | = |2 | = 1.

On en dduit que le shma est stable

Partie numrique.

on traitera le problme de diusion convection suivant :

u
u
2u
+c
2
t
x
x
u (x, 0)
u
(0, t)
x
u
(1, t)
x
avec (0

L2

= 0, x [0, 1] , t > 0
= u(0) (x)
= 0
= 0

< < 1)

0

(x)2 (1x)2
u(0) (x) = 16 sin 4 x
0.25
(0.25)4

La mise en uvre se droule ainsi :

1
on se xe N et on pose x =
, on xe t
N +1
On discrtise les conditions aux limites l'ordre 1 :

si

06x<

si

6x61

si

1<x61

uN +1 = uN

et

u0 = u1 .

On discrtise le problme de diusuion convection avec un schma de notre choix.

Par exemple le schma implicite dcentr amont

(n)
(3.10)
uj
(n)

(n1)
uj



 (n1)
 (n)
(n1)
(n)
(n1)
+ uj+1 , n > 0, 1 6 j 6 N

2uj
u uj1 +
u
2 j+1
x2 j1

(n)

uN +1 = uN

(n)

(n)

= u1

u0

Ce qui se met sous forme matricielle :

AU (n+1) = BU (n) avec A = et B =
3.9.

Exercices.

3.9.1.

Shma Lax-Wendro.

(3.11)

(n+1)
uj

 2 

 (n)
(n)
(n)
(n)
(n)

u uj1 +
uj+1 2uj + uj1
2 j+1
2

Erreur de troncature locale et globale

2
Stabilit L ,l
sous la condition CFL
Variation totale non conserve (voir sur dessin)

3.9.2.
(3.12)

(n)
uj

Shma Lax-Friedrichs.
(n+1)

uj

(n)
(n)

uj+1 + uj1  (n)
(n)

uj+1 uj1
2
2
2
et globale : noter que vj+1 + vj1 = 2vj + O (x )

Erreur de troncature locale

2
Stabilit L , stabilit l : condition CFL
Principe du maximum discret
Variation totale dcroissante

Universit de Pau et des Pays de l'Adour

17 Laboratoire de Mathmatiques Appliques

Janvier 2010

3.9.3.

Rsolution d'EDP par dirences nies

Schma de Newmark.

Pierre Puiseux

voir Quarteroni, p 473

4. Diffrences finies pour les problmes paraboliques

Introduction de la forme matricielle et des conditions aux limites

=R

+ condition

initiale
4.1.

Equation de la chaleur, pb modle.

u
2u
2
t
x
u (x, 0)
u (0, t)
u (1, t)

Maillage:

t > 0

et

= f

si

(x, t) ]0, 1[ ]0, +[

= u0 (x)
= 0
= 0

x > 0

Mise en quation du problme

condition initiale donne

sont le pas de temps et le pas d'espace, donns.

t,x = h = {(tn , xj ) = (nt, jx) , 1 j J

4.2.

et

n > 0}

Le shma d'Euler explicite.




u
t

n

n

u
x2

1
un+1
unj + O (t)
j
t




1
=
unj1 2unj + unj+1 + O x2
2
x
=

on rajoute les conditions aux limites de Dirichlet : ce qui fournit le schma :

un+1
= unj
j

t n
n
n

2u
+
u
+ [Link]
u
j1
j
j+1
2
x

unJ+1 = 0
un0 = 0
4.2.1.

Consistance.
(n)

(4.1)
4.2.2.

Stabilit L2 .

= O x2 t + O (t)

On oublie les conditions et on crit (4.1) sous la forme :

un+1 (x) = un (x) (un (x x) 2un (x) + un (x + x))

On prend la transforme de Fourier des deux membres :

(x) = un (x) 1 eix 2 + eix

un+1

d'o le coecient d'amplication :

1 eix 2 + eix
x
= 1 4 sin2
2

A () =

Schma
(4.2)

L2

stable si et seulement si pour tout

,|A ()| < 1

soit

t
1

2
x
2

Universit de Pau et des Pays de l'Adour

18 Laboratoire de Mathmatiques Appliques

Janvier 2010

4.2.3.

Rsolution d'EDP par dirences nies

Ecriture matricielle du shma. U

(n)

(n)
(n)
(n)
u1 , u2 , ..., uJ

Pierre Puiseux

t

U (n+1) = M U (n) + tF (n)

(4.3)

M =I

(4.4)

t
A
x2

et

..
..
.
.
.
1 2
.

.
.
.
..
..
..
A =
0
0
. .
.
..
..
..
2 1
0 0 1 2
= tridiagn (1, 2, 1)
Les valeurs propres de la matrice

(4.5)
donc les valeurs propres
(4.6)
4.2.4.

sont



k = 4 sin2 k x , 1 j J
2
de M sont


2
k = 1 4 sin k x , 1 j J
2

Stabilit, analyse matricielle.

Schma stable si et seulement si

(4.7)
On retrouve (4.2).
2

x 10

t 104

ou

k, |k | 1

soit

t
1

2
x
2

J = 100 n = 10000

pas de temps !

Sous la condition CFL, on a

(1) pour tout
(2)

k , k > 0

car

1 4 sin2 k 2 x > 1 2 sin2 k 2 x = sin (kx)

kM k = sup1kJ |k | = 1

soit encore

kM k = 1 4 sin

(4.8)


2

Pour l'inuence du second membre, voir exercice

4.2.5.

convergence.

sous la condition CFL, (

(n)

= erreur de troncature U (n) projection de

la solution exacte )

U (n+1) = M U (n) + tF (n) + (n)

(4.9)
en soustrayant

(4.9)

et

(4.3)

on obtient l'erreur

E (n+1)
E (n+1)

=
...
=

M E (n) + (n)
M n+1 E (0) +

M k (nk)

k=0,n1
Compte tenu de

E (0) = 0,

on en tire une majoration de l'erreur en norme euclidienne :

(n+1)
E

kM kk (nk)

k=0,n1

1
sup (n)
1 kM k n

Universit de Pau et des Pays de l'Adour

19 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Pierre Puiseux

(grace

2
t
t
x
2 t, et (n) = O (xt)+O x
= O (x) ,
1kM k = 4 x
2 sin
2
t
la condition CFL,
= O (1)), on obtient une majoration en O (x)
x2

4.2.6.

Convergence vers une solution stationnaire.

De plus

4.3.

Shma d'Euler implicite.

M U (n+1) = U (n) + tF (n)

(4.10)

M =I +

(4.11)
4.3.1.

Erreur de troncature locale et globale.

4.3.2.

Stabilit L2 , inconditionnelle.

4.3.3.

Stabilit l .

4.3.4.

Principe du maximum.

4.3.5.

Variation totale.

4.4.

TODO

Le teta-schma.

Le paramtre

t
A
x2

est pris dans l'intervalle

[0, 1],

ce schma consiste

mixer les deux schmas Euler implicite et explicite, et crire

n+1
un+1
unj
un+1
+ un+1
unj1 2unj + unj+1
j
j1 2uj
j+1
(4.12)

(1 )
= fjn+
2
2
t
x
x
n+
Dans cette criture on a pos fj
= f (jx, (n + ) t). Sous forme matricielle, tenant
compte des conditions aux limites de Dirichlet homognes :

U (n+1) = U n
4.4.1.

t
A U (n+1) + (1 ) U (n) + tF (n)
2
x

Consistance du schma.

on pose

tn+ = (n + ) t et vjn+ = v (jx, (n + ) t)

Proposition 18. Soit f une fonction dnies au voisinage d'un point a R, susament
rgulire, et soit [0, 1]. Alors

f (t + t) + (1 ) f (t t) = f (t) + (2 1) f 0 (t) t + O t2

Dmonstration.

DL de

f (t t)

Autrement dit, f (t + t) + (1 ) f (t
1
0
2 si =
(ou si f (t) = 0), d'ordre 1 sinon.
2

t)

est une approximation de

f (t),

d'ordre

Proposition 19. Le teta-schma est consistant d'ordre 1 en temps si 6= 12 et d'ordre 2

en temps si = 12 , et d'ordre 2 en espace. Pour =
Nicholson.
Dmonstration.

1
2

il s'appelle le schma de Cranck-

On crit l'erreur locale de troncature

n+1
= T1 + T2 + T3
j
avec

1
un+1
unj
j
t

un+1
un+1
+ un+1
=
j1 2
j
j+1
2
x


(1 ) n
n
n
=
u

2
u
+
u

j
j+1
j1
x2

T1 =
T2
T3

Universit de Pau et des Pays de l'Adour

20 Laboratoire de Mathmatiques Appliques

Janvier 2010

Rsolution d'EDP par dirences nies

Pierre Puiseux

Et on exprime


T1 =
puis en utilisant 18 la fonction

u
t

n+ 21

2u

x2

f () =

+ O t2

(x, )

n+1


2 u
+ O x2
=
2
x j
 2 n+ 21
 3 n+ 12 !




u
u
t
2
2
=
+
O
x
+
+
O
t
x2 j
2 x2 t j


T2

et enn

n


2 u
2
= (1 )
+
O
x
x2 j
 2 n+ 21
 3 n+ 12 !




u
u
t
= (1 )

+ O t2 + O x2
2
2
x j
2 x t j


T3

u tant la solution exacte, il vient :

 2 n+ 21
n+ 21
 3 n+ 12




u
u
t
u
=

+
(2 1)
+ O x2 + O t2
2
2
t j
x j
2
x t j




= (2 1) O (t) + O x2 + O t2

Finalement,

n+1
j


4.4.2.

Stabilit L2 .

en supposant

f = 0,

A () =

on trouve

1 4 (1 ) sin2 x
2
2 x
1 + 4 sin 2

|A ()| < 1 pour tout . Ce qui quivaut 1+4 (2 1) sin2 x

>
2
(s) = 1 + 4 (2 1) s > 0, s [0, 1].

et la stabilit a lieu lorsque

0, R ou encore
si > 21 cette ingalit est vrie
si < 21 cette ingalit est vrie si

et seulement si

(1) > 0

car

est dcroissante.

Ce qui a lieu pour



1
1
4

1
>
2

Le schma de Crank-Nicholson.

1
2
Consiste mixer les deux schmas implicite et explicite, et crire (ici sous forme

4.5.

C'est le

-schma

avec

matricielle)

(4.13)




un+1
unj
 n+1
j
n
n
n
= fjn
uj1 2un+1
+ un+1

j
j+1 + uj1 2uj + uj+1
2
t
2x

Sous forme matricielle, si l'on tient compte des conditions aux limites de Dirichlet
homognes :

(n+1)

= U A

U (n) + U (n+1)
2

+ tF (n)

soit

(n+1)

= I+

1

A
2

Universit de Pau et des Pays de l'Adour

A
2

1

U + I+ A
tF (n)
2
n

21 Laboratoire de Mathmatiques Appliques

Janvier 2010

4.6.

Rsolution d'EDP par dirences nies

Pierre Puiseux

1

tude matricielle de la stabilit. M = I + A

I

A
a pour valeurs propres
2
2


2k
2x
2

k = 2+k , avec k = 4 sin k 2 x Sp (A). La fonction (x) = 2+x

est dcrois
2
12 sin ( 2 x)
< 1...
sante car > 0 donc k est maximum pour 1 , et (M ) =
1+2 sin2 ( 2 x)
1
peauner car si >
les k ne sont pas toutes >0
2

Schma de Richardson (saute mouton).

(4.14)

un+1
un1

j
j
n
n
n
n

2u
+
u
j
j1 + fj
2t
x2 j+1

(4.15)

U (n+1) = U (n1) 2AU (n) + 2tF (n)

4.6.1.
4.6.2.

Consistance. Schma consistant d'ordre 2 en temps, et 2 en espace

t
Stabilit L2 . On pose = x
2
L2 -instable. Pour le montrer, on
(x, t) = (jx, nt) sous la forme

C'est un schma
schma au point

suppose que

f = 0,

et on crit le

u(n+1) (x) = u(n1) (x) + 2 u(n) (x + x) 2u(n) (x) + u(n) (x x)
on en prend la transforme de Fourier.



d
d
d
\
(n+1) () = u
(n1) () + 2 u
(n) () eix 2u
(n) () + u
(n) () eix
u\


d
(n1) () + 2 u
(n) () eix 2 + eix
= u\
x
d
(n1) () 8 u
(n) () sin2
= u\
2
qui se met sous forme matricielle :

(n+1) ()
u\
(n) ()
ud

!
= A ()

(n) ()
ud
(n1) ()
u\

avec

1
8 sin2 x
2
1
0

A () =
Le schma est

L2 -stable si et seulement si le spectre de A () est contenu dans le disque

|z| 6 1.
Les valeurs propres de

A ()

sont les racines du polynme



2 x
p () =
8 sin
() 1
2
x
= 2 + 8 sin2
1
2
q


2 x
2 x 2
qui admet les deux racines = 4 sin

1
+
4
sin

.
2
2
On a | | > 1 pour toute valeur de , ce qui montre l'instabilit inconditionnelle

du

schma.
4.7.

Le schma de Dufort-Franckel.

(4.16)
Consistance

C'est le schma Richardson stabilis

t n
n+1
n1
n
un+1
= un1
+
u

u
+
u
+
u
+ [Link]
j+1
j1
j
j
j
j
2
x
 2
t
:O
+ O (x2 ) (conditionnellement consistant).
x2

Stabilit : inconditionnellement stable en norme

Universit de Pau et des Pays de l'Adour

L2

22 Laboratoire de Mathmatiques Appliques