Discrétisation des problèmes paraboliques
Discrétisation des problèmes paraboliques
On a vu au paragraphe 1.3.2 comme exemple type de problème parabolique l’équation de la chaleur instationnaire :
ut − ∆u = f
qui fait intervenir la dérivée en temps d’ordre 1, ut , ainsi qu’un opérateur différentiel d’ordre 2 en espace. Pour
que ce problème soit bien posé, il faut spécifier des conditions aux limites sur la frontière de Ω, et une condition
initiale en t = 0.
Proposition 3.2 (Principe du maximum) Sous les hypothèses du théorème 3.1, soit u la solution du problème
(3.1.1) ;
1. si u0 (x) ≥ 0 pour tout x ∈ [0, 1], alors u(x, t) ≥ 0, pour tout t ≥ 0 pour tout x ∈]0, 1[.
2. 'u'L∞([0,1[×]0,T [) ≤ 'u'L∞ (]0,1[) .
67
3.2. DISCRÉTISATION PAR EULER EXPLICITE EN TEMPS. CHAPITRE 3. EDP PARABOLIQUES
Ces dernières propriétés peuvent être importantes dans le modèle physique ; supposons pas exemple que u repré-
sente une fraction massique. Par définition de la fraction massique, celle-ci est toujours comprise entre 0 et 1.
La proposition précédente nous dit que la quantité u donné par le modèle mathématique supposé représenter la
fraction massique d’une espèce qui diffuse dans un milieu, par exemple, est aussi comprise entre 0 et 1, dès que la
fraction massique initiale u0 est dans l’intervalle [0, 1] ce qui est plutôt une bonne nouvelle : le modèle mathéma-
tique respecte les bornes de la physique. Mais on ne peut pas en général calculer la solution de (3.1.1) de manière
analytique. On a recours à la disrétisation en temps et en espace pour se ramener à un système d’équations de
dimension finie. Il est souhaitable pour la validité du calcul que la solution approchée obtenue par la résolution de
ce système, qui est supposée approcher une fraction massique soit aussi comprise à tout instant entre 0 et 1. On dit
souvent d’une méthode de discrétisation (ou d’un schéma de discrétisation) qu’elle (ou il) est “robuste” ou “stable”
s’il préserve les bornes imposées par la physique (0 et 1 dans le cas de la fraction massique évoquée ci-dessus).
Pour calculer une solution approchée, on se donne une discrétisation en temps et en espace, qu’on notera D. On
choisit pour l’instant de discrétiser par différences finies en temps et en espace. La discrétisation consiste donc à se
donner un ensemble de points tn , n = 1, . . . , M de l’intervalle ]0, T [, et un ensemble de points xi , i = 1, . . . , N .
Pour simplifier, on considère un pas constant en temps et en espace. Soit : h = N1+1 = ∆x le pas de discrétisation
T
en espace, et k = ∆t = M , le pas de discrétisation en temps. On pose alors tn = nk pour n = 0, . . . , M et
xi = ih pour i = 0, . . . , N + 1. On cherche à calculer une solution approchée uD du problème (3.1.1) ; plus
précisement, on cherche à déterminer uD (xi , tn ) pour i = 1, . . . , N , et n = 1, . . . , M . Les inconnues discrètes
(n)
sont notées ui , i = 1, . . . , N et n = 1, . . . , M .
1 % (n) &
(n+1) (n)
(n) ūi − ūi (n) (n) (n)
R̃i = − ut (xi , tn ) et R̂i = 2ū i − ū i−1 − ū i+1 − uxx (xi , tn ).
k h2
Proposition 3.4 Le schéma (3.2.2) est consistant d’ordre 1 en temps et d’ordre 2 en espace, c’est à dire qu’il existe
C ∈ IR + ne dépendant que de u tel que :
(n)
|Ri | ≤ C(k + h2 ). (3.2.3)
' est (n)
Démonstration : On a vu lors de l’étude des problèmes elliptiques que l’erreur de consistance en espace R i
' (n)
d’ordre 2 (voir formule (2.2.30) page 18). Un développement de Taylor en temps donne facilement que Ri est
d’ordre 1 en temps.
3.2.2 Stabilité
On a vu à la proposition 3.2 page 67 que la solution exacte vérifie :
Si on choisit correctement les pas de temps et d’espcace, nous allons voir qu’on peut avoir l’équivalent discret sur
la solution approchée.
Définition 3.5 On dit qu’un schéma est L∞ -stable si la solution approchée est bornée dans L∞ indépendamment
du pas du maillage.
soit encore :
(n+1) (n) (n) (n)
ui = (1 − 2λ)ui + λui−1 + λui+1 .
3.2.3 Convergence
(n)
Définition 3.7 Soit u la solution du problème (3.1.1) et (ui ) i=1,...N, la solution de (3.2.2). On appelle erreur
n=1,...,M
de discrétisation au point (xi , tn ) la quantité eni = u(xi , tn ) − uni .
Théorème 3.8 Sous les hypothèses du théorème 3.1, et sous la condition de stabilité (3.2.4), il existe C ∈ IR + ne
dépendant que de u tel que
(n+1) (0)
'ei '∞ ≤ 'ei '∞ + T C(k + h2 ), pour tout i = 1, . . . , N et n = 0, . . . , M − 1.
(0) (n)
Ainsi, si 'ei '∞ = 0, alors maxi=1,...N 'ei ' tend vers 0 lorsque k et h tendent vers 0, pour tout n = 1, . . . M .
Le schéma (3.2.2) est donc convergent.
(n)
Démonstration : On note ūi = u(xi , tn ). On a donc, par définition de l’erreur de consistance,
(n+1) (n)
− ūi ūi 1 (n) (n) (n) (n)
− (2ūi − ūi−1 − ūi+1 ) = Ri . (3.2.5)
k h2
D’autre part, le schéma numérique s’écrit :
(n+1) (n)
− ui ui 1 (n) (n) (n)
− (2ui − ui−1 − ui+1 ) = 0. (3.2.6)
k h2
Retranchons (3.2.6) à (3.2.5), on obtient :
(n+1) (n)
ei − ei 1 (n) (n) (n) (n)
− (2ei − ei+1 − ei−1 ) = Ri ,
k h2
soit encore :
(n+1) (n) (n) (n) (n)
ei = (1 − 2λ)ei + λei−1 + λei+1 + kRi
(n) (n) (n)
Or (1 − 2λ)ei + λei−1 + λei+1 ≤ 'e(n) '∞ , car λ ≤ 21 , et donc comme le schéma est consistant, l’inégalité
(3.2.3) entraîne que :
(n+1)
|ei | ≤ 'e(n) '∞ + kC(k + h2 ).
On a donc par récurrence :
(n+1)
'ei '∞ ≤ 'e(0) '∞ + M kC(k + h2 )
ce qui démontre le théorème.
Donnons maintenant un exemple où lorsque la condition (3.2.3) n’est pas vérifiée, le schéma est instable.
u0 (x)
−1 −1 + ε 1
u (x) ≥ 0
0
u0 (x) )= 0 si x ∈] − 1; −1 + ε[
u0 (x) = 0 si x > −1 + ε
On considère le problème :
u − uxx = 0, ∀x ∈] − 1, 1[; ∀t > 0.
t
u(x, 0) = u0 (x), ∀x ∈] − 1, 1[ (3.2.7)
u(1, t) = u(−1, t) = 0, ∀t > 0.
On peut montrer que la solution exacte u de (3.2.7) vérifie u(x, t) > 0, ∀x ∈] − 1, 1[, ∀t > 0. En particulier,
(n)
pour un temps T > 0 donné, on a u(0, T ) > 0. Soit M ∈ IN et k = T /M . Soit ui la solution approchée par
(3.2.2), sensée approcher u(xi , tn ) (i ∈ {−N, . . . , N }, n ∈ N ). On va montrer que uM0 = 0 pour k et h choisis
de manière non admissible ; ceci montre que le schéma ne peut pas converger. Calculons uM 0 :
M−1
uM
0 = (1 − 2λ)u0 + λuM−1
−1 + λuM−1
1 .
Donc uM
0 dépend de
uM
0 )→ u(0, T ).
u − unum = u − uD + uD − unum .
On sait que l’erreur de discrétisation u − uD tend vers 0 lorsque h et k tendent vers 0, sous condition de stabilité
(3.2.4), c.à.d.
1
λ≤ .
2
Pour contrôler l’erreur entre la solution uD du schéma (3.2.2) et la solution numérique obtenue unum , on cherche
à estimer l’amplification de l’erreur commise sur la donnée initiale. Rappelons que le schéma s’écrit :
(n+1) (n) (n) (n)
ui = (1 − 2λ)ui + λui−1 + λui+1 ,
k
avec λ = . Ce schéma se met sous la forme u(n+1) = AU (n) , avec
h2
1 − 2λ λ 0 ... 0
. . ..
λ 1 − 2λ λ . .
. .. . .. . ..
A= 0 0
.. .
. .. λ 1 − 2λ λ
0 ... 0 λ 1 − 2λ
Définition 3.9 (Stabilité au sens des erreurs d’arrondi) Supposons que l’on commette une erreur ε0 sur la condi-
'0 , s’écrit donc u
tion initiale. La nouvelle condition initiale u '0 = u0 + ε0 . A cette nouvelle condition initiale cor-
(n) (n) (n)
respond une nouvelle solution calculée u ' = u + ε . On dit que le schéma est stable au sens des erreurs
d’arrondi s’il existe C > 0 indépendant de n tel que ε(n) ≤ Cε(0) .
On peut trouver une condition suffisante pour que le schéma 3.2.2 soit stable au sens des erreurs d’arrondi. En
effet, on va démontrer le résultat suivant :
k
Proposition 3.10 On suppose que λ = h2 < 21 . Alors le schéma 3.2.2 est stable au sens des erreurs d’arrondi.
Démonstration : Soit donc une condition initiale perturbée u '0 = u0 + ε0 à laquelle on associe une nouvelle
(n) (n) (n) (n) n 0
solution calculée u ' = u + ε . On a ε = A ε . Comme A est symétrique, A est diagonalisable dans
IR. Soient µ1 , . . . , µN les valeurs propres de A, et e1 , . . . , eN les vecteurs propres associés, c’est–à–dire tels que
Aei = µi ei , ∀i = 1, . . . N . On décompose la perturbation ε0 sur la base des vecteurs propres :
N
. N
.
ε0 = ai ei . On a donc An ε0 = ai µni ei = ε(n) .
i=1 i=1
sup |µi | ≤ 1
i=1...N
Soit VP(A) l’ensemble de valeurs propres de A. Alors VP(A) = {1 + λµ, µ ∈ VP(B))}. Or VP(B) =
{−4 sin2 2(Njπ+1) , j = 1, . . . , N } (voir Lemme 3.11 plus loin). Pour que ε(n) < ε0 , il faut donc que :
jπ
sup |1 − 4λ sin2 | < 1,
j=1,...,N 2(N + 1)
c.à.d.
jπ 1
λ sin2 < .
2(N + 1) 2
Une condition suffisante pour avoir une diminution de l’erreur est donc que λ < 12 .
Lemme 3.11 (Valeurs propres de B) L’ensemble VP(B) des valeurs propres de la matrice B définie par (3.2.8)
est donné par :
jπ
VP(B) = {−4 sin2 , j = 1, . . . , N }.
2(N + 1)
Démonstration : Les valeurs propres B peuvent se calculer à partir des valeurs propres de l’opérateur continu ; on
commence donc par chercher u solution de :
/
−u$$ + αu = 0,
u(0) = u(1) = 0.
et donc
(Bv (k) )i = sin kπ(i − 1)h − 2 sin kπih + sin kπ(i + 1)h
Problème continu avec conditions aux limites périodiques On considère le problème avec conditions aux
limites périodiques
ut − u(xx) = 0, t ∈]0, T [, x ∈]0, 2π[,
u(0, t) = u(2π, t), ∀t ∈]0, T [, (3.2.9)
u(x, 0) = u0 (x).
Le problème (3.2.9) est bien posé, au sens où ∀u0 ∈ C([0, 2π]), il existe une unique u ∈ C 2 (]0, 2π[×]0, T [, IR)
solution de (3.2.9). On suppose que u0 ∈ L2 (]0, 2π[). On rappelle que L2 est un espace de Hilbert, et que
{einx , n ∈ ZZ } est une [Link] 3 de L2 (]0, 2π[). On décompose donc la condition initiale dans cette
base hilbertienne : u0 (x) = cn (0)einx (au sens de la convergence dans L2 ). Dans un premier temps, calcu-
n∈ZZ
lons formellement les solutions de (3.2.9) sous la forme d’un développement dans la base hilbertienne :
.
u(x, t) = cn (t)einx .
n∈ZZ
ut − uxx = 0.
La condition de périodicité est bien vérifiée par u donnée par (3.2.10). Enfin on a bien : u(x, t) → u0 (t) lorsque
t → 0, donc la condition initiale est vérifiée. On peut remarquer qu’il y a “amortissement” des coefficients de
Fourier cn (0) lorsque t augmente, c.à.d. qu’on a : cn (t) ≤ cn (0), ∀ t > 0.
On prend comme condition initiale u0 (x) = ap eipx , pour p ∈ ZZ fixé . En discrétisant, on obtient : u0j (x) =
ap eipjh , pour j = 1, . . . , N, avec h = N2π 0 0
+1 . On a bien u0 = uN +1 = 0. Calculons :
(1)
uj = (1 − 2λ)ap eipjh + λap eip(j−1)h + λap eip(j+1)h
(1)
donc : uj = ap eipjh ξp . On appelle ξp le facteur d’amplification associé à la fonction eipx (appelé aussi “p-ième
mode”). On a donc :
(1) (0)
uj = ξp uj
..
.
u(n) = (ξ )n u(0)
j p j
Calculons ξp :
ξp = 1 − 2λ + 2λ cos ph
ph
= 1 − 2λ + 2λ(1 − 2 sin2 )
% & 2
p
= 1 − 4λ sin2 N2π
+1 , 2 .
2 3
2 2π p 1
Pour avoir |ξp | < 1, il faut λ sin , < Une condition suffisante pour que le schéma soit stable au
N +1 2 4
1
sens de Von Neumann est que : λ < . Remarquons que c’est la même condition que pour la stabilité des erreurs
2
d’arrondis.
(n)
. 2 . 2
(u − uD )j ≤ cp (0)(e−p nk
− ξpn )eipjh + cp (0)(e−p nk
− ξpn )eipjh
|p|≤A |p|≥A
On a donc : . . 2
(n)
(u − uD )j ≤X +2 |cp (0)|, avec X = |cp (0)|(e−p nk
− ξpn )
|p|≥A |p|≤A
.
et 2 |cp (0)| ≤ 2ε. Montrons maintenant que X → 0 lorsque h → 0. Remarquons que
|p|≥A
2 2 ph
e−p nk
− ξpn = e−p T
− ξpn , et ξp = 1 − 4λ sin2 .
2
ph p 2 h2 k ph
Or, sin2 = + O(h4 ), et λ = 2 . Donc : 4λ sin2 = p2 k + O(kh2 ). On en déduit :
2 4 h 2
2 3T /k 2 3
ph T ph
(ξp )n = 1 − 4λ sin2 et donc ln ξpn = ln 1 − 4λ sin2 = −T p2 + O(h2 ).
2 k 2
2
On en déduit que ξpn → e−p T
lorsque h → 0. Tous les termes de X tendent vers 0, et X est une somme finie ;
(n)
on a donc montré que (u − uD )j tend vers 0 lorsque h tend vers 0.
Remarque 3.13 On peut adapter la technique de Von Neumann au cas Dirichlet homogène sur [0, 1], en effec-
tuant le développement de u par rapport aux fonctions propres de l’opérateur u$$ avec conditions aux limites de
Dirichlet : .
u(x, t) = cn (t) sin(nπx).
L’avantage du développement en série de Fourier est qu’il marche pour n’importe quel opérateur linéaire à condi-
tion d’avoir pris des conditions aux limites périodiques.
Soit k un pas (constant) de discrétisation , on rappelle que les schémas d’Euler explicite et implicite pour la
discrétisatision de ce problème s’ecrivent respectivement :
y (n+1) − y (n)
Euler explicite : = f (y (n) ), n ≥ 0 (3.3.13)
k
y (n+1) − y (n)
Euler implicite : = f (y (n+1) ), n ≥ 0, (3.3.14)
k
avec y (n) = y0 . On rappelle également que le θ-schéma, où θ est un paramètre de l’intervalle [0, 1] sécrit :
Remarquons que pour θ = 0 on retrouve le schéma (3.3.13) et pour θ = 1 le schéma (3.3.14). On peut facilement
adapter le θ schéma à la résolution des équations paraboliques. Par exemple, le θ-schéma pour la discrétisation en
temps du problème (3.1.1), avec une discrétisation par différences finies en espace s’écrit :
(n+1) (n)
ui −ui (n+1) (n+1) (n+1) (n) (n) (n)
k = hθ2 (−2ui + ui−1 + ui+1 ) + 1−θ
h2 (−2ui + ui−1 + ui+1 ), ; n ≥ 0, i = 1 (3.3.16)
(0)
ui = u0 (xi ), i = 1, . . . , N.
1
Si θ = 0, on retrouve le schéma d’Euler explicite ; si θ = 1, celui d’Euler implicite. Dans ce cas où θ = ce
2
4 5
schéma s’appelle schéma de Crank -Nicolson . Notons que dès que θ > 0, le schéma est implicite, au sens où on
(n+1) (n)
n’a pas d’expression explicite de ui en fonction des uj .
(n+1) (n)
(n) ūj − ūj θ % (n+1) (n+1) (n+1)
& 1−θ %
(n) (n) (n)
&
Rj = + −2ū i + ū i−1 + ū i+1 + −2u i + u i−1 + u i+1
k h2 h2
4. John Crank (6 February 1916 –3 October 2006) mathématicien britannique.
5. Phyllis Nicolson (21 Septembre 1917 – 6 Octobre 1968) mathématicienne britannique.
(n,1) k
Tj = (ūt )(xj , tn ) + (utt )(xj , tn ) + R1 avec |R1 | ≤ Ck 2 .
2
(n,2)
Faisons de même pour Tj :
(n,2) 0 1
Tj = θ ūxx (xj , tn+1 ) + R2 avec |R2 | ≤ Ch2 .
Or ūxx (xj , tn+1 ) = ūxx (xj , tn ) + kūxxt (xj , tn ) + R3 avec |R3 | ≤ Ck 2 , donc :
(n,2) 0 1
Tj = θ uxx (xj , tn ) + +kuxxt(xj , tn ) + R4 avec |R4 | ≤ C(h2 + k 2 ).
(n,3)
De même pour Tj , on a :
(n,3)
Tj = (1 − θ)uxx (xj , tn ) + R5 , avec |R5 | ≤ Ck 2 .
(n) k ∂
Rj = ut (xj , tn ) − uxx (xj , tn ) ut (xj , tn ) + θk(uxx )(xj , tn ) + R
2 ∂t
avec R = R1 + R4 + R5
1
• Si θ = , on a un schéma d’ordre 2 en temps et en espace.
2 0 61 17
En effet, k2 (ūtt )(xj , tn ) − θk(ūxxt )(xj , tn ) = ∂t
∂
k 2 (ūt )(xj , tn ) − θ(ūxx )(xj , tn ) et ut − uxx = 0.
1
• Si θ )= : on a un schéma d’ordre 2 en espace et d’ordre 1 en temps.
2
1
Proposition 3.15 (Stabilité au sens de Von Neumann) Si θ ≥ le θ-schéma est inconditionnellement stable.
2
1
En particulier, les schémas d’Euler implicite et de Crank-Nicolson sont inconditionnellement stables. Si θ < le
2
schéma est stable sous condition.
1
λ≤ .
2(1 − 2θ)
(On retrouve en particulier que le schéma d’Euler explicite n’est que si λ ≤ 21 ).
(0)
En discrétisant la condition intiale (mode de Fourier) on obtient uj = eipjh et on cherche le facteur d’amplifcation
ξp tel que u1j = ξp u0j = ξp eipjh ; en appliquant le schéma ci–dessus pour n = 0, on obtient :
et donc :
1 − 2λ(1 − θ) + 2λ(1 − θ) cos ph 1 − 4λ(1 − θ) sin2 ph/2
ξp = =
(1 + 2λθ) − 2λ cos ph 1 + 4λθ sin2 ph
2
Pour que le schéma soit stable au sens de Von Neumann, il faut que : |ξp | < 1 pour tout p, soit encore :
ph ph
1 − 4λ(1 − θ) sin2 < 1 + 4λθ sin2 (3.3.18)
2 2
et
ph ph
− 1 < 1 + 4λθ sin2
4λ(1 − θ) sin2 (3.3.19)
2 2
L’inégalité (3.3.18) est toujours vérifiée. En ce qui concerne l’inégalité (3.3.19), on distingue deux cas :
1
1. Si θ ≤ 2 alors 0 ≤ 1 − θ ≤ θ et dans ce cas (3.3.19) est toujours vraie.
1
2. Si θ < 2 , on veut : 9 :
2 ph 2 ph
4λ (1 − θ) sin − θ sin <2
2 2
Il faut donc que
; <−1
1 2 ph
λ< (1 − 2θ) sin
2 2
Une condition suffisante est donc :
1 1
λ≤ si θ < .
2(1 − 2θ) 2
On rappelle que ce schéma est inconditionnellement stable au sens de Von Neumann. On va montrer de plus qu’il
est L∞ –stable :
(n)
Proposition 3.16 (Stabilité L∞ pour Euler implicite) Si (uj )j=1,...,N est solution du schéma (3.3.20), alors :
de même :
(n+1) (n) (0)
min uj ≥ min uj ≥ min uj (3.3.22)
j=1,...,N j=1,...,N j=1,...,N
(n+1) (n)
On en déduit : uj0 ≤ maxj=1,...,N uj , ce qui prouve que
(n+1) (n)
max uj ≤ max uj .
j=1,...,N j=1,...,N
Alors 'e(n+1) '∞ ≤ 'e(0) '∞ + T C(k + h2 ). Si 'e(0) '∞ = 0, le schéma est donc convergent (d’ordre 1 en temps
et 2 en espace.
et donc :
'e(n+1) '∞ ≤ 'eh '∞ + kC(k + h2 )
On en déduit, par récurrence sur n, que :
Si le domaine est rectangulaire, ce problème se discrétise facilement à l’aide de θ schéma en temps et de différences
finies en espace, en prenant un maillage rectangulaire. On peut montrer, comme dans le cas 1D, la consistance, la
stabilité, la L∞ stabilité, la stabilité au sens de Von Neumann
Montrer que u est bien définie de [0, 1] × IR $+ dans IR et est solution de (3.5.23) au sens suivant :
u ∈ C ∞ ([0, 1] × IR $+ , IR),
ut (x, t) − uxx (x, t) = 0, ∀x ∈ [0, 1], ∀t ∈ IR $+ ,
(3.5.25)
u(0, t) = u(1, t) = 0, ∀t ∈ IR $+ ,
'u(., t) − u0 'L2 (]0,1[) → 0, quand t → 0.
2. Montrer qu’il existe une unique fonction u solution de (3.5.25).
Exercice 24 (Exemple de schéma non convergent) Suggestions en page 92, corrigé en page 95
Soit u0 ∈ L2 (] − 4, 4[). On note u l’unique solution (au sens vu en cours ou en un sens inspiré de l’exercice
précédent) du problème suivant :
1. Soient T > 0, v ∈ C 1 ([0, 1], IR + ), a0 , a1 ∈ IR et u0 ∈ C([0, 1]). On considère le problème d’évolution suivant :
ut (x, t) − uxx (x, t) + v(x)ux (x, t) = 0, x ∈]0, 1[, t ∈]0, T [,
u(0) = a0 , u(1) = a1 , (3.5.28)
u(x, 0) = u0 (x).
dont on cherche à approcher la solution par différences finies. On choisit pour cela le schéma de la question 1 de
l’exercice 5 pour la discrétisation en espace, et on discrétise par le schéma d’Euler implicite en temps avec un pas
de temps uniforme k = PT où P ≥ 1.
(p) (p)
1.1 Ecrire le schéma ainsi obtenu et montrer qu’il admet une solution qu’on notera U = (ui ) i=1,...N, , où ui est
p=1,...,P
censé être une approximation de u(xi , tp ), où tp = pk, p = 0, . . . , P.
1.2 Montrer que
(p)
min(min u0 , min(a0 , a1 )) ≤ ui ≤ max(max u0 , max(a0 , a1 )), pour tout i = 1, . . . , N et p = 1, . . . , P .
[0,1] [0,1]
dont on cherche à approcher la solution par différences finies. On choisit pour cela le schéma de la question 2 de
l’exercice 6 pour la discrétisation en espace, et on discrétise par le schéma d’Euler implicite en temps avec un pas
de temps uniforme k = PT où P ≥ 1.
(p) (p)
2.1 Ecrire le schéma ainsi obtenu et montrer qu’il admet une solution qu’on notera U = (ui ) i=1,...N, , où ui est
p=1,...,P
censé être une approximation de u(xi , tp ), où tp = pk, p = 0, . . . , P.
(p)
2.2 Montrer que si a0 ≥ 0, a1 ≥ 0 et u0 ≥ 0, alors on a : ui ≥ 0, pour tout i = 1, . . . , N et p = 1, . . . , P .
3. On considère maintenant Ω =]0, 1[2 ; soient v ∈ C ∞ (Ω, (IR + )2 ) (v(x) est donc un vecteur de IR 2 ), a ∈
C(∂Ω, IR) et u0 ∈ C(Ω, IR + ). En s’inspirant des schémas étudiés aux questions 3 et 4, donner une discrétisation
où u0 et ε sont donnés (ε > 0). On reprend dans la suite de l’exercice les notations du cours.
1. Donner un schéma d’approximation de (3.5.34) différences finies à pas constant et centré en espace et Euler
explicite à pas constant en temps. Montrer que l’erreur de consistance est majorée par C(k+h2 ) ; avec C dépendant
de la solution exacte de (3.5.34). Sous quelle(s) condition(s) sur k et h a-t’on 'un '∞ ≤ 'u0 '∞ , ∀n ≤ M ?
Donner un résultat de convergence pour ce schéma.
2. Même question que 1. en remplaçant Euler explicite par Euler implicite.
3. En s’inspirant du schéma de Crank-Nicolson (vu en cours) construire un schéma d’ordre 2 (espace et temps).
Sous quelle(s) condition(s) sur k et h a-t’on 'un '2 ≤ 'u0 '2 , ∀n ≤ M ? Donner un résultat de convergence pour
ce schéma.
4. Dans les schémas trouvés aux questions 1., 2. et 3. on remplace l’approximation de ux par une approximation
décentrée à gauche. Quel est l’ordre des schémas obtenus et sous quelle(s) condition(s) sur k et h a-t’on 'un '∞ ≤
'u0 '∞ où 'u'u0'2 , ∀n ≤ M ? Donner un résultat de convergence pour ces schémas.
Soit u0 une fonction donnée de [0, 1] dans IR. On s’intéresse ici à la discrétisation du problème suivant :
un+1
0 = un+1 0
N +1 = 0, n ∈ IN ; ui = u0 (xi ), = i ∈ {0, . . . , N + 1}. (3.5.39)
Pour n ∈ N , on note un = (un1 , . . . , unN )t ∈ IR N .
1. (Consistance) Soit T > 0. Pour n ∈ IN, et i ∈ {1, . . . , N }, on note Rin l’erreur de consistance (définie en
cours) du schéma numérique (3.5.37), (3.5.39) [resp. du schéma numérique (3.5.38), (3.5.39)]. Montrer qu’il
existe C ∈ IR, ne dépendant que de u et T , t. q. |Rin | ≤ C(k + h2 ), pour tout i ∈ {1, . . . , N } et tout n ∈ IN,
t.q. kn ≤ T .
Ecrire une discrétisation espace-temps du problème (3.5.26) de l’exercice 23 avec le schéma de Crank-Nicolson
en temps, et par volumes finis avec un maillage uniforme en espace, en utilisant un schéma décentré amont pour le
terme d’ordre 1 ux (x, t).
(schéma de Dufort-Frankel).
k
3. Montrer que le schéma (3.5.42) est consistant avec (3.5.40) quand h, k → 0 sous la condition h → 0.
On suppose que u0 ∈ C(]0, 1[, IR) . On rappelle que dans ce cas, il existe une unique fonction u ∈ C 2 (]0, 1[
×]0, T [, IR) ∩ C([0, 1] × [0, T ], IR) qui vérifie (3.5.43). On cherche une approximation de la solution de ce
problème, par une discrétisation par différences finies en espace et en temps. On se donne un ensemble de points
{tn , n = 1, . . . , M } de l’intervalle ]0, T [, et un ensemble de points {xi , i = 1, . . . , N }. Pour simplifier, on
considère un pas constant en temps et en espace. Soit : h = N1+1 le pas de discrétisation en espace, et k = M T
, le
pas de discrétisation en temps. On pose alors tn = nk pour n = 0, . . . , M et xi = ih pour i = 0, . . . , N + 1. On
cherche à calculer une solution approchée uapp du problème (3.5.43) ; plus précisement, on cherche à déterminer
(n)
uapp (xi , tn ) pour i = 1, . . . , N , et n = 1, . . . , M . Les inconnues discrètes sont notées ui , i = 1, . . . , N et
n = 1, . . . , M .
On considère le schéma suivant :
;
1 % (n+1) (n) (n−1)
& 1 (n+1) (n+1) (n+1) i = 1, . . . N,
3u i − 4u i + u i + (2u i − u i−1 − u i+1 ) = 0,
2k h 2 n = 1, . . . , M,
0
ui = u0 (xi ), i = 1, . . . , N, (3.5.44)
u1 = u (x ), i = 1, . . . , N,
i 1 i
(n) (n)
u0 = uN +1 = 0, ∀n = 1, . . . , M,
U n+1 = BW n (3.5.45)
2 −1 0 ··· ··· 0
n+1 −1 2 −1 0 ··· 0
u1
2k 0 −1 2 −1 ··· 0
où U n+1 = ... , B = (3Id + 2 A)−1 , A = . .. .. .. .. ,
.. (3.5.46)
h .. . . . .
.
n+1
uN
0 ··· 0 −1 2 −1
0 ··· ··· 0 −1 2
n−1 n
u1 u1
n n−1 .. n ..
et W ne dépend que de U = . et U = . . (3.5.47)
un−1
N unN
On se propose, dans cet exercice, de montrer l’existence d’une solution faible au problème (3.5.48)-(3.5.50), à
partir de l’existence de la solution approchée donnée par un schéma numérique. L’inconnue de ce problème est la
fonction u de [0, 1] × [0, T ] dans IR, elle doit être solution des équations suivantes :
∂u ∂ 2 ϕ(u)
(x, t) − (x, t) = v(x, t), x ∈]0, 1[, t ∈]0, T [, (3.5.48)
∂t ∂x2
∂ϕ(u) ∂ϕ(u)
(0, t) = (1, t) = 0, t ∈]0, T [, (3.5.49)
∂x ∂x
u(x, 0) = u0 (x), x ∈]0, 1[, (3.5.50)
> >
∂ψ ∂2ψ
(u(x, t) (x, t) + ϕ(u(x, t)) 2 (x, t) + v(x, t)ψ(x, t))dxdt + u0 (x)ψ(x, 0)dx = 0,
D ∂t ∂x ]0,1[
(3.5.52)
∀ψ ∈ CT2,1 (IR 2 ),
1. Soit a > 0, Pour s ∈ IR, on pose ga (s) = s + aϕ(s). Montrer que ga est une application strictement
croissante bijective de IR dans IR.
2. Soit w = (w i )i=1,...,N ∈ IR N . On pose w 0 = w 1 etw N +1 = wN . Montrer qu’il existe un et un seul couple
(u, w) ∈ IR N × IR N , u = (ui )i=1,...,N , w = (wi )i=1,...,N , t.q. :
On peut donc définir une application F de IR N dans IR N par w /→ F (w) = w où w est solution de
(3.5.55)–(3.5.56).
3. On munit IR N de la norme usuelle ' · '∞ . Montrer que l’application F est strictement contractante. [On
pourra utiliser la monotonie de ϕ et remarquer que, si a = ϕ(α) et b = ϕ(β), on a |α − β| ≥ (1/L)|a − b|,
où L ne dépend que de ϕ.]
5. Soit w = (wi )i=1,...,N t.q. w = F (w). Montrer que pour tout i ∈ {1, . . . , N } il existe un+1
i ∈ IR t.q.
wi = ϕ(un+1
i . Montrer que {un+1
i , i = 1, . . . , N } est solution de (3.5.54).
et
N
. +1
(ϕ(un+1
i ) − ϕ(un+1 2
i+1 )) ≤ C2 h, pour tout n ∈ {0, . . . , M }. (3.5.59)
i=0
Dans la suite de l’exercice, il s’agit de passer à la limite (quand N, M → ∞) pour trouver une solution de (3.5.48)-
(3.5.50).
√
Pour M ∈ IN$ donné, on prend N = M 2 (et donc h et k sont donnés et k = T h), on définit (avec les uni trouvés
dans les questions précédentes) une fonction, uh , sur [0, 1] × [0, T ] en posant
Question 6 Montrer que les suites (uh )M∈IN ! et (ϕ(uh ))M∈IN ! sont bornées dans L∞ (]0, 1[×]0, T [) (on rappelle
que h est donné par M ).
pour tout η ∈ IR $+ , avec ϕ(uh )(·, t) prolongée par 0 hors de [0, 1].
2. 'ϕ(uh )(·, t) − ϕ(uh )(·, s)'L2 (]0,1[) ≤ C|t − s|, pour tout t, s ∈ [0, T ].
Etudier ensuite les racines de léquation r2 − αr − 1 = 0 et montrer que l’une de ses racines est, en module,
supérieure à 1.
4. Reprendre la méthode développée à la question 2, en montrant que l’équation caractéristique pour ξ est mainte-
nant :
p(r) = ar2 + br + c = 0,
avec
1 1 2 cos(ph) 1 1
+ 2, b = −
a= 2
et c = 2 − .
2k h h h 2h
Etudier ensuite les racines de cette équation.
2 2 2
Comme (x, t) → e−n π t an sin(nπt) est continue (pour tout n ∈ IN∗ ), on en déduit que u est continue sur
IR×]ε, ∞[, et finalement sur IR×]0, ∞[ car ε > 0 est arbitraire.
Pour dériver terme à terme la série définissant u, il suffit également d’obtenir sur ]ε, ∞[×IR (pour tout ε > 0) une
majoration du terme général de la série des dérivées par le terme général d’une série convergente (indépendant de
(x, t) ∈ IR×]ε, ∞[. On obtient cette majoration en remarquant que, pour (x, t) ∈ IR×]ε, ∞[,
2
π 2 t2 2
π 2 ε2
| − n2 π 2 e−n an sin(nπx)| ≤ n2 π 2 e−n r2 'u0 '2
et ceci donne ut = uxx sur IR × IR ∗t et donc aussi un [0, 1] × IR ∗+ . Le fait que u(0, t) = u(1, t) pour tout t > 0 est
immédiat car sin nπt = sin 0 = 0, pour tout n ∈ IN ∗ .
Il reste à montrer que u(., t) → u0 dans L2 (]0, 1[) quand t → 0.
avec
N
. 2
π 2 t2
u(N ) (x, t) = an e−n sin(nπx)
n=1
N
.
(N )
u0 (x) = an sin(nπx).
n=1
(N )
Il est clair que, pour tout N ∈ IN∗ , on a u(N ) (., t) → u0 uniformément sur IR, quand N → ∞, et donc
(N )
u(N ) (., t) → u0 dans L2 (]0, 1[).
Comme
∞
. .∞
1 2 2 2 1 (N )
'u(., t) − u(N ) (., t)'22 = a2n e−2n π t ≤ a2n = 'u0 − u0 '22 → 0
2 2
n=N +1 n=N +1
2
quand N → ∞, on en déduit que u(., t) → u0 , qdt → 0, dans L (]0, 1[).
2) On note w la différence de 2 solutions de (3.5.25). On a donc
w ∈ C ∞ ([0, 1] × IR ∗+ , IR)
wt − wxx = 0 sur [0, 1] × IR ∗+
w(0, t) = w(1, t) = 0 pour t > 0
w(., t) → 0, dans L2 (]0, 1[), quand t → 0
Soit 0 < ε < T < ∞. On intègre l’équation wwt − wwxx = 0 sur ]0, 1[×]ε, T [. En utilisant une intégration par
parties (noter que w ∈ C ∞ ([0, 1] × [ε, T ]), on obtient :
> 1 > 1 > 1 > T
1 2 1 2
w (x, T )dx − w (x, ε).dx + wx2 (x, t)dxdt = 0.
2 0 2 0 0 ε
D’où l’on déduit 'w(., T )'2 ≤ 'w(., ε)'2 . Quand ε → 0, on a 'w(., ε)'2 → 0, on a donc 'w(., T )'2 = 0 et donc,
comme w(., t) est contenue sur [0, 1], wX ∈ [0, 1]. Comme T > 0 est arbitraire, on a finalement
N +1 N +1
u01 = u0 (ih, 0), i=− ,...,
2 2
Soit maintenant n ∈ {0, . . . , M }. On a :
N +1 N +1
un+1
i =0 pour i=− et i=
2 2
k 0 n 1 N +1 N +1
un+1
i 2
= uni +
ui+1 + uni−1 − 2uni , i = − + 1, . . . , − 1.
h 2 2
2) On va montrer, par récurrence (finie) sur n, que :
; <
N +1 N +1 N +1
Si n ∈ {0, . . . , M + 1}, i ∈ − ,..., et i ≥ − + n alors uni = 0. (3.7.60)
2 2 4
N +1
Pour initialiser la récurrence, on suppose que n = 0 et i ≥ − . On a alors
4
N +1 8
ih ≥ − = −2 > −3
4 N +1
et donc u0i = 0.
Soit maintenant n ∈ {0, . . . , M }. On suppose que?l’hypothèse de récurrence
@ est vérifiée jusqu’au rang n, et on
démontre la propriété au rang n + 1. Soit donc i ∈ − N2+1 , . . . , N2+1 tel que i ≥ − N4+1 + (n + 1). Alors :
– Si i = N2+1 on a bien uNi
+1
= 0.
– Si i < 2 , les indices i − 1, i et i + 1 sont tous supérieurs ou égaux à − N4+1 + n, et donc par hypothèse de
N +1
récurrence, 2 3
n+1 2k k k
ui = ui 1 − 2 + 2 uni+1 + 2 uni−1 = 0.
n
h h h
1 8
On a donc bien démontré (3.7.60). On utilise maintenant l’hypothèse k = h, c’est–à–dire M+1 = N +1 . On a alors
N +1
− + M + 1 = −2(M + 1) + M + 1 = −(M + 1) < 0.
4
On en déduit que si n ∈ {0, . . . , M + 1} et i ≥ 0, alors i ≥ − N4+1 + n. On en déduit que uni = 0 pour
? @
n ∈ {0, . . . , M + 1} et i ∈ 0, . . . , N2+1 . On remarque alors que
; ; << ; ; <<
N +1 N +1 N +1
max |uM+1 i − ū M+1
i |, i ∈ − , . . . , ≥!la max |ū M+1
i |, i ∈ 0, . . . ,
2 2 2
k2
ūn+1
i = ūni + kut (ih, nk) + utt (ξ1 , t1 ), (3.7.62)
2
h2
ūni−1 = ūni − hux (ih, nk) + uxx (ξ2 , t2 ), (3.7.63)
2
h2 h3 h4
ūni−1 = ūni − hux (ih, nk) + uxx (ih, nk) − uxxx (ih, nk) − uxxxx(ξ3 , t3 ), (3.7.64)
2 6 24
h2 h3 h4
ūni+1 = ūni + hux (ih, nk) + uxx (ih, nk) + uxxx (ih, nk) + uxxxx(ξ4 , t4 ). (3.7.65)
2 6 24
On en déduit :
k h
Rin = ut (ih, nk) + utt (ξ1 , t1 ) + αux (ih, nk) + α uxx (ξ2 , t2 )
2 2
h2
−µuxx (ih, nk) − µ (uxxxx(ξ3 , t3 ) + µuxxxx(ξ4 , t4 )) ,
24
et donc, comme u est solution de (3.5.32), pour h assez petit, on a :
|Rin | ≤ C1 (h + k),
où C1 ne dépend que de u. Le schéma (3.5.33) est donc consistant d’ordre 1 en temps et en espace.
(b) Cherchons les conditions pour que un+1i s’écrive comme combinaison convexe de uni , uni−1 et uni+1 .
On peut réécrire le schéma (3.5.33) :
αk 2µk µk αk µk
un+1
i = auni + buni+1 + cuni−1 , avec a = 1 − − 2 , b = 2 et c = + 2.
h h h h h
h2
k≤ . (3.7.66)
αh + 2µ
Si h et k vérifient la condition (3.7.66), on pose : M n = maxi=1...N uni (resp. mn = mini=1...N uni .
Comme un+1 i est une combinaison convexe de uni , uni−1 et uni+1 , on a alors : un+1
i ≤ M n ∀i =
n+1 n n+1 n n+1
1, . . . , N (resp. ui ≥m ∀i = 1, . . . , N ) et donc : M ≤ M (resp. m ≥ mn ). On a ainsi
montré que :
'un+1 '∞ ≤ 'un '∞ .
On a de même :
'un '∞ ≤ 'un−1 '∞ .
..
.
Donc, sous la condition (3.7.66), on a 'un+1 '∞ ≤ 'un '∞ et donc 'un '∞ ≤ 'u0 '∞ , pour tout
n = 1, . . . , N.
(c) En retranchant l’égalité (3.7.61) au schéma (3.5.33), on obtient l’équation suivante sur eni :
1 n+1 α µ
(ei − eni ) + (eni − eni−1 ) − 2 (eni−1 − 2eni + eni+1 ) = Rin .
k h h
ce qu’on peut encore écrire :
kα kµ kµ
en+1
i = (1 − − 2 2 )eni + eni−1 2 + kRin . (3.7.67)
h h h
Sous la condition de stabilité (3.7.66), on obtient donc :
|en+1
i | ≤ 'en+1 '∞ + C1 (k + h)k,
n
|ei | ≤ 'en−1 '∞ + C1 (k + h)k,
.. .. ..
. ≤ . + .
|e1i | ≤ 'e0 '∞ + C1 (k + h)k,
Si à t = 0, on a 'e0 ' = 0, alors on éduit des inégalités précédentes que |eni | ≤ C1T (k + h) pour tout
n ∈ IN. Le schéma est donc convergent d’ordre 1.
2. Schéma explicite centré.
h2 h3
ūni−1 = ūni − hux (ih, nk) + uxx(ih, nk) − uxxx(ξ5 , t5 ),
2 6
h2 h3
ūni+1 = ūni + hux (ih, nk) + uxx (ih, nk) + uxxx (ξ6 , t6 ),
2 6
on obtient maintenant :
k h2
Rin = ut (ih, nk) + utt (ξ1 , t1 ) + αux (ih, nk) + α (uxxx(ξ5 , t5 ) + µuxxx (ξ6 , t6 ))
2 12
h2
−µuxx (ih, nk) − µ (uxxxx(ξ3 , t3 ) + µuxxxx(ξ4 , t4 )) ,
24
On en déduit que
|Rin | ≤ C3 (k + h2 ),
où C3 = max( 12 'utt '∞ , 16 'uxxx'∞ , 12
1
'uxxxx'∞ ).
(b) Le schéma s’écrit maintenant :
2µk µk αk µk αk
un+1
i = ãuni + b̃uni+1 + c̃uni−1 , avec ã = 1 − , b̃ = 2 − et c̃ = 2 + .
h2 h h h h
Remarquons que l’on a bien : ã + b̃ + c̃ = 1. Pour que un+1 i soit combinaison convexe de uni , uni+1 et
uni−1 , il faut et il suffit donc que ã ≥ 0, b̃ ≥ 0, et c̃ ≤ 0. L’inégalité c̃ ≥ 0 est toujours vérifiée. Les
deux conditions qui doivent être vérifiées par h et k s’écrivent donc :
2µk
i. ã ≥ 0, i.e. 1 − h2 ≥ 0, soit encore
h2
k≤ .
2µ
µk αk
ii. b̃ ≥ 0 i.e. h2 − h ≥ 0, soit encore
µ
h≤ .
2α
Le schéma centré est donc stable sous les deux conditions suivantes :
µ 1 2
h≤ et k ≤ h . (3.7.68)
2α 2µ
Pour obtenir une borne d’erreur, on procède comme pour le schéma (3.5.33) : on soustrait la définition
de l’erreur de consistance au schéma numérique, et on obtient :
en+1
i = ãeni + b̃eni+1 + c̃eni−1 + kRin . (3.7.69)
Par le même raisonnement que pour le schéma décentré, on obtient donc que si e0i = 0, on a |eni | ≤
C4 (k + h2 ), avec C4 = T C3 .
h2 h3 h4
u(x + h, t) = u(x, t) + hux (x, t) + uxx (x, t) + u(3) (x, t) + u(4) (α, t),
2 6 24
h2 h3 $$$ h4 (4)
u(x − h, t) = u(x, t) − hux (x, t) + uxx (x, t) − u (x, t) + u (β, t),
2 6 24
k2
u(x, t + k) = u(x, t) + kut (x, t) + utt (x, τk ), τk ∈ [t, t + k].
2
De ces développements de Taylor, il ressort que l’erreur de consistance vérifie |R| ≤ C(k + h2 ), où C ne dépend
que de u. Le schéma est donc explicite d’ordre 1 en temps et 2 en espace.
Cherchons alors les conditions pour que :
'un '∞ ≤ 'u0 '∞ .
Par définition,
'un '∞ = max |unj |.
j=1,...,N
n+1 n
On essaye d’abord de vérifier que : 'u '∞ ≤ 'u '∞ , c’est-à-dire :
max |un+1
j | ≤ max |unj |,
j=1,...,N j=1,...,N
min un+1
j ≥ min unj .
j=1,...,N j=1,...,N
mn+1 ≥ mn
..
.
où unj est calculé par (3.7.70), et ūnj est la solution du problème (3.5.34) en xj = jh et tn = nk.
On a donc, par définition de l’erreur de consistance,
1 n+1 1 n ε
(ū − ūnj ) + (ū − ūnj−1 ) − 2 (−2ūnj + ūnj+1 + ūnj−1 ) = Rjn
k j 2h j+1 h
où |Rin | ≤ C(k + h2 )
ce qui entraine :
1 n+1 1 n ε
(e − enj ) + (e − enj−1 ) − 2 (−2enj + enj+1 + enj−1 ) = Rjn
k j 2h j h
soit encore : 2 3 2 3 2 3
n+1 2εk n k kε j+1 εk k
ej = 1 − 2 ej + − + e + + en + kRjn .
h 2h h2 h2 2h j−1
Analyse numérique II, Télé-enseignement, M1 100 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
de même que précédemment, on obtient sous les conditions (3.7.71)
|en+1
j | ≤ 'en '∞ + C(k + h2 )k
..
.
|enj | ≤ 'en−1 '∞ + C(k + h2 )k
..
.
|e1j | ≤ 'e0 '∞ + C(k + h2 )k.
Et donc sous les conditions (3.7.71) on a 'en '∞ qui tend vers 0 lorsque k, h → 0, ce qui prouve que le schéma est
convergent.
Exercice 28 page 85
(n)
1. Notons Ri l’erreur de consistance en (xi , tn ). Pour le schéma (3.5.37), on a donc par définition :
(n+1) (n)
(n) ū −ū (n+1) (n+1) (n+1)
Ri = i k i + 1
h2 (2ūi − ūi−1 − ūi+1 ) − ūn+1
i
(n)
= R̃i + R̂in ,
où
ūn+1 − ūni
R̃in = i
− ut (xi , tn ) est l’erreur de consistance en temps
k
et
1
(2ūn+1
R̂in =
i − ūn+1 n+1
i−1 − ūi+1 ) − (uxx (xi , tn )) est l’erreur de consistance en espace.
h2
On a vu (voir (2.2.30)) que 5 4 5
5 5 h2 5∂ u 5
5 n5 5
5R̂i 5 ≤ sup (·, tn )55 , ∀i ∈ {1, . . . , N }
12 [0,1] 5 ∂x4
Effectuons maintenant un développement de Taylor en fonction du temps d’ordre 2 :
k2
u(xi , tn+1 ) = u(xi , tn ) + kut + utt (xi , ξn )
2
u(xi , tn+1 ) − u(xi , tn ) k
avec ξn ∈ [tn , tn+1 ]. Donc − ut = utt (xi , ξn ). Comme ξn ∈ [0, T ], et utt admet un
k 2
maximum (à xi fixé) dans [0, T ] (qui est compact), on a donc
5 5 k
5 n5
5R̃i 5 ≤ max |utt (xi , ·)| .
2 [0,T ]
Analyse numérique II, Télé-enseignement, M1 101 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
Par conséquent,
5 5 5 5 5 5 k 5 4 5
5 n n5 5 n 5 5 n5 h2 5∂ u 5
n
|Ri | = 5R̃i + R̂i 5 ≤ 5R̃i 5 + 5R̂i 5 ≤ max |utt (xi , ·)| + max 5 4 (·, tn+1 )55 .
5
2 [0,T ] 12 [0,1] ∂x
Donc |Rin | ≤ C(k + h2 ) avec
2 3
1 1 ∂ 4u
C= max 'utt 'L∞ ([0,1]×[0,T ] ; ' 4 'L∞ ([0,1]×[0,T ] .
2 6 ∂x
Le calcul de l’erreur de consistance pour le schéma (3.5.38) s’effectue de manière semblable.
2. Le schéma (3.5.37) est complètement implicite alors que le schéma (3.5.38) ne l’est que partiellement, puisque le
terme de réaction est pris à l’instant n. Le schéma (3.5.37) s’écrit : AU n+1 = U n avec U n+1 = (U1n+1 , . . . UN
n+1 t
) , Un =
n n t
(U1 , . . . , UN ) , et
1 + 2λ − k −λ 0 ... 0
.
−λ 1 + 2λ − k −λ . . 0
0
A=
..
. 0
0 0 −λ 1 + 2λ − k
k
où λ = 2 . Notons que par définition, A est symétrique. De même, le schéma (3.5.38) s’écrit : BU n+1 = U n
h
avec
1 + 2λ −1 0 ... 0
..
−1 1 + 2λ . 0
1
0 −1
B=
1+k
0 0 −1 1 + 2λ
On a donc A = λAh , où Ah est définie en (2.2.26) page 15, avec ci = 1−k λ
λ , et B = k+1 Ah avec ci = λ1 . Dans
les deux cas, les matrices sont donc s.d.p. en vertu de la proposition 2.2 page 15. Notons que l’hypothèse k ∈]0, 1[
est nécessaire dans le cas du premier schéma, pour assurer la positivité de ci .
3. Le schéma (3.5.38) s’écrit
(1 + k)uni = un+1
i + λ(2un+1
i − un+1 n+1
i+1 − ui−1 ).
On montre facilement par récurrence que maxj=1,...,N unj ≤ (1 + k)n maxj=1,...,N u0j , (voir preuve de la stabilité
(n) (0)
L∞ d’Euler implicite page 80) et que min uj ≥ (1 + k)n min uj . On en déduit que
j=1,...,N j=1,...,N
Analyse numérique II, Télé-enseignement, M1 102 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
On en déduit le résultat, avec C1 (T ) = eT . De même, pour le schéma (3.5.37), on montre par récurrence que :
1
'u(n) '∞ ≤ 'u(0) '.
(1 − k)n
Mais pour k ∈]0, α[, avec α ∈]0, 1[, on a :
1 1
≤ 1 + βk, = avec β = .
(1 − k) (1 − α)
On en déduit par un calcul similaire au précédent que
(1 − k)T /k ≤ eβT ,
1
'e(n+1) ' ≤ 'e(n) ' + kC(k + h2 ).
1−k
Par récurrence sur n, on obtient alors
2 3n
1
'e(n) '∞ ≤ [kC(k + h2 ) + 'e0 '∞ ]
1−k
d’où
'e(n) '∞ ≤ C2 (T, α)(T C(k + h2 ) + 'e0 '∞ ).
De même, pour le schéma (3.5.38), on écrit l’erreur de consistance :
(n+1) (n+1) (n+1)
ūj − ūnj uj+1 + ūn+1
j−1 − 2ūj (n,2)
− − ūnj = Rj
k h2
et donc :
(n+1) (n+1) (n+1) (n) (n,2)
ej (1 + 2λ) − λej−1 − λej+1 = ej (1 + k) + kRj .
Par des raisonnements similaires à ceux de la question 3 on obtient alors :
(n)
'ej ' ≤ (1 + k)n ('e(0) ' + kC(k + h2 ))
d’où
'e(n) '∞ ≤ C1 (T )('e(0) ' + kC(k + h2 )).
Analyse numérique II, Télé-enseignement, M1 103 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
Exercice 30 page 86
1. On s’intéresse ici à l’ordre du schéma au sens des différences finies. On suppose que u ∈ C 4 ([0, 1] × [0, T ]) est
solution de (3.5.40) et on pose
ūn+1
j − ūjn−1 ūnj−1 − 2ūnj + ūnj+1
Rjn = − , j = 1, . . . , N − 1, k = 1, . . . , M − 1.
2k h2
On cherche une majoration de Rjn en utilisant des développement de Taylor. Soit j ∈ {1, . . . , N − 1}, k ∈
{1, . . . , M − 1}. Il existe (ξi , ti ) ∈ [0, 1] × [0, T ], i = 1, . . . , 4, t.q. :
k2 k3
ūn+1
j = ūnj + kut (jh, nk) + utt (jh, nk) + uttt (ξ1 , t1 ),
2 6
k2 k3
ūn−1
j = ūnj − kut (jh, nk) + utt (jh, nk) − uttt (ξ2 , t2 ),
2 6
h2 h3 h4
ūnj−1 = ūnj − hux (jh, nk) + uxx (jh, nk) − uxxx (jh, nk) − uxxxx(ξ3 , t3 ),
2 6 24
h2 h3 h4
ūnj+1 = ūnj + hux (jh, nk) + uxx (jh, nk) + uxxx(jh, nk) + uxxxx(ξ4 , t4 ).
2 6 24
On en déduit :
k2 h2
Rjn = ut (jh, nk) + (uttt (ξ1 , t1 ) + uttt (ξ2 , t2 )) − uxx (jh, nk) − (uxxxx(ξ3 , t3 ) + uxxxx(ξ4 , t4 )) ,
12 24
et donc, comme u est solution de (3.5.40),
|Rjn | ≤ C1 (k 2 + h2 ),
Analyse numérique II, Télé-enseignement, M1 104 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
2
Le facteur d’amplification est donc, pour tout t ∈ IR + , le nombre e−p t . Ce facteur est toujours, en module,
inférieur à 1. On va maintenant chercher la solution du schéma numérique sous la forme :
où ξ0 et ξ1 ∈ IR sont donnés (ils donnent u0j et u1j pour tout j ∈ ZZ ) et ξn ∈ IR est à déterminer de manière à ce
que la première équation de (3.5.41) ) soit satisfaite.
Ce facteur ξn va dépendre de k, h et p. Pour k et h donnés, le schéma est stable au sens de Von Neumann si, pour
tout p ∈ IR, la suite (ξn )n∈IN est bornée. Dans le cas contraire, le schéma est (pour ces valeurs de k et h) dit
instable au sens de Von Neumann.
Un calcul immédiat donne que la famille des unj , définie par (3.7.72), est solution de la première équation si et
seulement si la suite (ξn )n∈IN vérifie (on rappelle que ξ0 et ξ1 sont donnés) :
ξn+1 − ξn−1 2
= 2 (cos ph − 1)ξn , n ≥ 2,
2k h
ou encore, en posant
4k
α= (cos ph − 1) (≤ 0),
h2
ūn+1
j − un−1
j ūnj−1 − (ūn+1
j + ūn−1
j ) + ūj+1
Sjn = − 2
, j = 1, . . . , N − 1, k = 0, . . . , M − 1.
2k h
En reprenant la technique de la question 1, il existe (ξi , ti ), i = 1, . . . , 6 t.q.
h2 h2 k2 k2
Sjn = (uttt (ξ1 , t1 ) + uttt (ξ2 , t2 )) − (uxxxx (ξ3 , t3 ) − uxxxx (ξ4 , t4 )) + 2 htt (ξ5 , t5 ) + 2 utt (ξ6 , t6 ).
12 24 2h 2h
Ce qui donne, avec C2 ne dépendant que de u,
2 3
k2
|Sjn | ≤ C2 h2 + k 2 + 2 , j = 1, . . . , N − 1, k = 0, . . . , M − 1.
h
k
Le schéma est donc consistant quand h → 0 avec → 0.
h
Analyse numérique II, Télé-enseignement, M1 105 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
4. On reprend la méthode développée à la question 2, la suite (ξn )n doit maintenant vérifier la relation suivante
(avec ξ0 , ξ1 donnés).
ξn+1 − ξn−1 2 cos(ph) ξn−1 + ξn+1
= ξn − ,n≥2
2k h2 h2
c’est à dire :
2 3 2 3
1 1 2 cos(ph) 1 1
ξn+1 + − ξn + ξn−1 − = 0, n ≥ 2.
2k h2 h2 h2 2k
L’équation caractéristique est maintenant :
p(r) = ar2 + br + c = 0,
avec
1 1 2 cos(ph) 1 1
+ ,b=− a= et c = 2 − .
2k h2 h2 h 2h
Pour montrer la stabilité au sens de Von Neumann, il suffit d’après (3.7.74) de montrer que les deux racines du
polynôme p sont de module inférieur ou égal à 1. On note r1 et r2 ces deux racines (qui peuvent être confondues)
et on distingue 2 cas :
1. 1er cas : Les racines de p ne sont pas réelles. Dans ce cas, on a |r1 | = |r2 | = γ et
c
γ = | | < 1,
a
car k > 0.
2. 2ème cas : Les racines de p sont réelles. Dans ce cas, on remarque que
c
r1 r2 = < 1,
a
2 2 cos ph
et l’une des racines, au moins, est donc entre −1 et 1 (strictement). De plus on a p(1) = 2 − ≥0
h h2
2 2 cos ph
et p(−1) = 2 + ≥ 0, l’autre racine est donc aussi entre -1 et 1 (au sens large).
h h2
On en déduit que le schéma (3.5.42) est stable au sens de Von Neumann.
∂u ∂ 2 ϕ(u)
(x, t) − (x, t) = v(x, t), x ∈]0, 1[, t ∈]0, T [, (3.7.76)
∂t ∂x2
∂ϕ(u) ∂ϕ(u)
(0, t) = (1, t) = 0, t ∈]0, T [, (3.7.77)
∂x ∂x
u(x, 0) = u0 (x), x ∈]0, 1[, (3.7.78)
Analyse numérique II, Télé-enseignement, M1 106 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
1. T > 0, v ∈ L∞ (]0, 1[×]0, T [),
2. ϕ croissante, lipschitzienne de IR dans IR,
3. u0 ∈ L∞ (]0, 1[) et ϕ(u0 ) lipschitzienne de [0, 1] dans IR.
Un exemple important est donné par ϕ(s) = α1 s si s ≤ 0, ϕ(s) = 0 si 0 ≤ s ≤ L et ϕ(s) = α2 (s − L) si s ≥ L,
avec α1 , α2 et L donnés dans IR $+ . Noter pour cet exemple que ϕ$ = 0 sur ]0, L[.
Les ensembles ]0, 1[ et D =]0, 1[×]0, T [ sont munis de leur tribu borélienne et de la mesure de Lebesgue sur cette
tribu.
On appelle “solution faible" de (3.5.48)-(3.5.50) une solution de :
> >
∂ψ ∂2ψ
(u(x, t) (x, t) + ϕ(u(x, t)) 2 (x, t) + v(x, t)ψ(x, t))dxdt + u0 (x)ψ(x, 0)dx = 0,
D ∂t ∂x ]0,1[
(3.7.80)
∀ψ ∈ CT2,1 (IR 2 ),
où ψ ∈ CT2,1 (IR 2 ) signifie que ψ est une fonction de IR 2 dans IR deux fois continûment dérivable par rapport à x,
∂ψ ∂ψ
une fois continûment dérivable par rapport à t et t.q. (0, t) = (1, t) = 0, pour tout t ∈ [0, T ] et ψ(x, T ) = 0
∂x ∂x
pour tout x ∈ [0, 1]}.
Analyse numérique II, Télé-enseignement, M1 107 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
et comme u vérifie (3.5.49), on a
∂ϕ(u) ∂ϕ(u)
(0, t) = (1, t) = 0, t ∈]0, T [.
∂x ∂x
En tenant compte de ces relations et en ré–intégrant par parties, on obtient :
> >
∂ 2 ϕ(u) ∂ 2 ψ(u)
2
(x, t)ψ(x, t)dxdt = − ϕ(u)(x, t) (x, t)dxdt. (3.7.83)
D ∂x D ∂x2
Analyse numérique II, Télé-enseignement, M1 108 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
3. Soit w 1 et w2 ∈ IR N et soit w1 = F (w 1 ) et w2 = F (w 2 ). Par définition de F , on a :
2k 1 k0 1
u1i − u2i + 2
(wi − wi2 ) = 2 (w 1i−1 + w 1i+1 ) − (w2i−1 + w2i+1 ) , pour tout i = 1, . . . , N. (3.7.84)
h h
Comme ϕ est monotone, le signe de wi1 − wi2 = ϕ(u1i ) − ϕ(u2i ) est le même que celui de u1i − u2i , et donc
2k 1 2k
|u1i − u2i + (w − wi2 )| = |u1i − u2i | + 2 |wi1 − wi2 |. (3.7.85)
h2 i h
Et comme ϕ est lipschitzienne de rapport L, on a
d’où :
1 1
|u1i − u2i | ≥ |w − wi2 |. (3.7.86)
L i
On déduit donc de (3.7.84),(3.7.85) et(3.7.86) que
1 1 2k k0 1
|wi − wi2 | + 2 |wi1 − wi2 | ≤ 2 |w 1i−1 − w 2i−1 | + |w 1i+1 − w 2i+1 | , pour tout i = 1, . . . , N.
L h h
On a donc
1
|wi1 − wi2 | ≤ h2
max |w 1i − w2i |, pour tout i = 1, . . . , N.
1 + 2kL i=1,...,N
1
d’où on déduit que 'w1 − w2 '∞ ≤ C'w1 − w 2 '∞ avec C = h2
< 1. L’application F est donc bien
1+ 2kL
strictement contractante.
4. Soit {un+1 i , Si {un+1
i , i = 1, . . . , N } est solution de (3.5.54) et wi = ϕ(un+1 i ) pour i ∈ {1, . . . , N },
alors on remarque que (un+1 i )i=1,...,N et (w )
i i=1,...,N vérifient (3.5.55)–(3.5.56) avec w i = wi pour i =
1, . . . , N. On en déduit que w = F (w).
5. Soit w = (wi )i=1,...,N t.q. w = F (w). Montrer que pour tout
Par définition de F , on a F (w) = w̃ avec (ũ, w̃) ∈ IR N × IR N , ũ = (ũi )i=1,...,N , w̃ = (w̃i )i=1,...,N , t.q. :
Analyse numérique II, Télé-enseignement, M1 109 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
Supposons que i est tel que un+1
i = minj=1,...,N j u n+1
. Comme ϕ est croissante, on a dans ce cas : ϕ(un+1
i−1 ) −
n+1 n+1 n+1 n+1 n
ϕ(ui ) ≥ 0 et ϕ(ui+1 − ϕ(ui ) ≥ 0, et on en déduit que minj=1,...,N uj ≥ ui − kB d’où, par hypothèse
de récurrence, minj=1,...,N un+1
j ≥ −A − nkB − kB. Un raisonnement similaire en considérant maintenant i
n+1 n+1 n+1 n
tel que ui = maxj=1,...,N uj conduit à : maxj=1,...,N uj ≤ ui + kB ≤ A + nkB + kB. On a donc
bien :−A − (n + 1)kB ≤ un+1 i ≤ A + (n + 1)kB, pour tout i = 1, . . . , N et tout n = 0, . . . , M .
On en déduit alors que 'un 'L∞ (]0,1[) ≤ cu0 ,v,T , avec cu0 ,v,T = A + BT.
a2 b2
En utilisant l’inégalité a2 − ab = 2 − 2 , on obtient :
N
1 . n 2
An ≥ αn+1 − αn , avec αn = (u ) .
2k i=1 i
En développant Bn , on obtient :
N N
1 0. n+1 n+1 n+1
. 1
Bn = − 2 (ϕ(ui−1 ) − ϕ(ui ))ui + (−ϕ(un+1
i ) + ϕ(un+1 n+1
i+1 ))ui ) .
h i=1 i=1
N −1
1 0. 1
Bn = 2
(ϕ(un+1 n+1
i+1 ) − ϕ(ui ))(un+1 n+1
i+1 − ui ) .
h i=1
Enfin, on majore Cn :
Bcu0 ,v,T
Cn ≤ .
h
L’égalité An + Bn = Cn entraîne donc :
N −1
1 . Bcu0 ,v,T
αn+1 − αn + (ϕ(un+1 n+1 2
i+1 ) − ϕ(ui )) ≤ .
Lh2 i=1 h
Analyse numérique II, Télé-enseignement, M1 110 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
En sommant pour n = 0 à M − 1, et en notant que αM ≥ 0,on obtient alors :
M−1 N −1
1 . . Bcu0 ,v,T
(ϕ(un+1 n+1 2
i+1 ) − ϕ(ui )) ≤ + α0 .
Lh2 n=0 i=1 h
h 2
Il reste à remarquer que α0 ≤ 2k cu0 ,v,T pour conclure que :
M−1
. N.−1
h 1
(ϕ(un+1 n+1 2
i+1 ) − ϕ(ui )) ≤ C1 , avec C1 = Lcu0 ,v,T (B + cu0 ,v,T ).
n=0 i=1
k 2
An + Bn = Cn , (3.7.88)
N N
. un+1 − un . ϕ(un+1 n+1
i−1 ) − 2ϕ(ui ) + ϕ(un+1
i+1 )
avec An = i i
(ϕ(un+1
i ) − ϕ(uni )), Bn = − (ϕ(un+1
i )−
i=1
k i=1
h2
N
.
ϕ(uni )) et Cn = vin (ϕ(un+1
i ) − ϕ(uni )).
i=1
En utilisant le caractère lipschitzien de ϕ, on obtient la minoration suivante :
N −1
1 .
An ≥ (ϕ(un+1
i ) − ϕ(uni ))2 . (3.7.89)
Lk i=1
En développant Bn , on obtient :
N N
1 0. n+1 n+1 n+1
. 1
Bn = − 2
(ϕ(u i−1 ) − ϕ(u i ))(ϕ(u i ) − ϕ(u n
i )) + (−ϕ(un+1
i ) + ϕ(un+1 n+1
i+1 ))(ϕ(ui ) − ϕ(uni ))) .
h i=1 i=1
1 %. &
N −1 .N
n+1 n+1 n+1
Bn = − 2 n
(ϕ(ui ) − ϕ(ui+1 ))(ϕ(ui+1 ) − ϕ(ui+1 )) + (−ϕ(un+1
i ) + ϕ(un+1 n+1
i+1 ))(ϕ(ui ) − ϕ(uni ) .
h i=0 i=1
N −1
1 0. 1
Bn = 2
(ϕ(un+1 n+1
i+1 ) − ϕ(ui ))((ϕ(un+1 n+1
i+1 ) − ϕ(ui )) − (ϕ(uni+1 ) − ϕ(uni )) .
h i=1
a2 b2
En utilisant à nouveau la relation a(a − b) ≥ 2 − 2 , on obtient :
N −1
1 .
Bn ≥ βn+1 − βn , avec βn = (ϕ(uni+1 ) − ϕ(uni ))2 (3.7.90)
2h2 i=1
Analyse numérique II, Télé-enseignement, M1 111 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010
3.7. CORRIGÉS DES EXERCICES CHAPITRE 3. EDP PARABOLIQUES
Enfin, on majore Cn par :
N N
1 . . k
Cn ≤ (ϕ(un+1
i ) − ϕ(uni ))2 + C k≥C . (3.7.91)
2Lk i=1 i=1
h
d’autre part, en utilisant que le fait que le premier terme est positif, on obtient par (3.7.92) une majoration sur βM ,
et donc sur βn pour tout n ≤ M :
C
βn ≤ + β0 . (3.7.94)
h
Il ne reste donc plus qu’à majorer β0 pour obtenir (3.5.58) et (3.5.59). Par définition, on a
N
. −1
ϕ(u0i ) − ϕ(u0i+1 )
β0 = .
i=1
2h2
En utilisant le fait que ϕ est lipschitzienne et que la différence entre u0i et u0i+1 est en h, on obtient (3.5.59) à partir
de (3.7.93) et (3.5.58) à partir de (3.7.94).
ce qui prouve que la suite (uh )M∈IN ! est bornée dans L∞ (]0, 1[×]0, T [). Comme ϕ est continue, on en déduit
immédiatement que (ϕ(uh ))M∈IN ! est bornée dans L∞ (]0, 1[×]0, T [)
Analyse numérique II, Télé-enseignement, M1 112 Université Aix-Marseille 1, R. Herbin, 27 septembre 2010