Processus stochastiques et équations aux dérivés partielles

CHAPITRE II :

Généralités sur les EDPs et la

méthode des différences finie

0. Introduction:

Le métier d'ingénieur est toujours besoin d'utilisation des logiciels de modélisation. Ces logiciels résolvent des équations telles que les équations aux dérivées partielles, mais la résolution ici est faite par une méthode discrète. Donc ces équations différentielles ne peuvent en général pas être résolues de façon exacte. Elles sont résolues de façon approchée, à l'aide des méthodes numériques.

Les méthodes numériques ne donnent pas la solution véritable du problème que l'on cherche à résoudre. Des méthodes numériques mal employées peuvent conduire à des résultats totalement faux, Il est indispensable pour un ingénieur de posséder des notions de base sur les méthodes numériques utilisé pour résoudre son problème.

Dans ce chapitre on présente une méthode numérique pour la résolution des EDPs,

La méthode des différences finies.

1. Généralité sur les équations aux dérivées partielles :

1.1. Définitions :

· Equation aux dérivées partielles :

Une équation aux dérivées partielles fait intervenir plusieurs variables indépen- dantes (temps et espace pour les équations de l'ingénieur), ainsi que les dérivées partielles de la variable dépendante ( c.-à-d. la solution recherchée) par rapport à ces variables indépendantes.

Exemple : l'équation de convection (parfois appelée advection)

=

ax 0 (2.1)

ac

+
at +u

ac

est une EDP, la variable dépendent est c , les variables indépendantes sont le temps t et l'espace x. La grandeur u (homogène à une vitesse) peut (ou non)être fonction de t ,x et c.

· Ordre d'une EDP :

L'ordre d'une EDP est défini exactement de la même façon que pour une EDO : c'est l'ordre le plus élevé parmi toutes les dérivées partielles de l'EDP.

Exemple : l'EDP (2.1) est une EDP d'ordre 1 (par rapport à t ou x ). En revanche, l'EDP suivant :

ac

+
at +u

ac

+ D

ax

ax2 = 0 (2.2)

a²c

est une EDP d'ordre 2 car sa dérivée partielle d'ordre le plus élevé est une dérivée seconde(en l'occurrence par rapport à x).

· EDP linéaires, quasi-linéaires et non- linéaires :

EDP linéaire : est une EDP qui ne fait intervenir que des combinaisons linéaires des dérivées partiales de la variable dépendante.

Exemple : l'équation suivante est linéaire :

a^2u A+2B ax^2

a^2u ~~~ ~~~~ ~ ~ ~2.3~ ~~~ ~ ~ ~~, ~, ~~ ~~ , ~~ ay) = 0

où les coefficients A, B et C sont dépend de X et Y (connues), F est une fonction lineaire et U est la fonction inconnue (à déterminer).

EDP quasi-linéaire : est une équation linéaire par rapport aux dérivées partielles d'ordres supérieurs.

Exemple : l'équation précédente (2.3) est quasi-linéaire si les coefficients A, B et C sont en fonction de U et/ou de X, Y et F peut être non linéaire.

EDP non-linéaire : est une EDP où l'une des dérivées partielles intervient comme argument d'une fonction non-linéaire.

Exemple : l'équation suivante:

Co

~~~~ ~~~~ ~ ~~ ~~~~ ~2.4~ ~~~~ ~ ~~ ~~~ ~~~ ~ ~~ ~~~ ~~~ ~ 0

est non-linéaire si au moins des fonctions Co(x) et C1(x) n'est pas linéaire par rapport à x.

1.2. Conditions initiales et conditions aux limites :

Contrairement aux EDOs, les conditions initiales ne suffisent pas à assurer l'unicité de solution. Il faut également fournir des conditions aux limites. Pour certains types d'équations (ex. EDPs elliptiques), le concept de condition initiale n'a pas de sens.

Les conditions initiales et les conditions aux limites se distinguent de la manière suivante :

· Une condition initiale s'applique pour une valeur donnée (et unique) d'une variable indépendante. A partir de cette condition initiale, il est possible de déduire la solution pour toutes les autres valeurs de la variable indépendante.

· Une condition aux limites est appliquée en tout point de la frontière du domaine sur lequel on souhaite résoudre l'équation (et non en un point unique). Il n'est pas possible de déterminer la solution en partant simplement d'un seul point de la limite de domaine et en progressant à l'intérieur de celui-ci, car la solution est également conditionnée par sa valeur en tous les autres points de la frontière.

1.3. Classification des EDPs du second ordre :

Les EDPs du second ordre représente une classe importante des EDPs du monde de l'ingénierie, on traite dans quel suit des EDPs second ordre quasi-linéaires de la forme générale :

a²u A + B
aX²

a²u a²u

+

aY2 D

au

+ E

aX

aY + F = 0 (2.5)

au

où A ,B,C,D,F sont des fonctions de X,Y et u. X et Y sont les variables
indépendantes (ce pourrait être t,x,y etc) de l'EDP et u la variable dépendante.

Selon la valeur des coefficients A ,B et C, l'EDP est dite hyperbolique, parabolique ou elliptique.

EDPs hyperbolique :

Une EDP du type (2.5) est dite hyperbolique si son discriminant A= B² -- 4AC est strictement positif.

Exemple : l'équation suivante est du type hyperbolique :

a²u

at²

A. ax2

² = 0 (2.6)

En effet, l'équation (2.5) peut être mise sous la forme (2.6) en posant X = t, Y = x, A = 1, C = --A.², et B = D = E = F = 0. Il est facile de vérifier que A= B² -- 4AC > 0. A noter que A. est homogène à une vitesse.

EDPs paraboliques :

Une EDP du type (2.5) est dite paraboliques si son discriminant A= B² -- 4AC est nul.

Exemple : l'équation suivante est du type paraboliques :

ax2 = 0, v > 0 (2.7)

au

-- v

at

a²u

En effet, l'équation (2.5) peut être mise sous la forme (2.7) en posant X = t, Y = x , D = 1, C = v, et A = B = E = F = 0. On vérifiera que A= B² -- 4AC = 0.

EDPs elliptique :

Une EDP du type (2.5) est dite elliptique si son discriminant A= B² -- 4AC est strictement négatif.

Exemple : l'équation suivante est du type elliptique :

0y2 = Q (2.8)

a²u ax2 +

a²u

En effet, l'équation (2.5) peut être mise sous la forme (2.8) en posant X = t, Y = x , A = C = 1, et B = D = E = F = 0, on peut vérifier que A= B² -- 4AC < 0.

· Une aide mnémotechnique :

L'<< astuce >> suivante peut etre utilisé pour déterminer la nature d'une EDP du seconde ordre : il suffit, dans l'équation originale, de remplacer les dérivées a^Pu / a XP (p = 1,2) par X^P et a^Pu / a Y^P (p = 1,2) par Y^Pet le second membre par une constante. Ainsi (2.5) devient :

AX²+BXY + Y² + DX + EX + F = Cst (2.9)

L'équation (2.9) est l'équation d'une courbe conique dans le plan(X, Y). Si cette courbe est une ellipse, l'EDP est elliptique ; si la courbe est une parabole, l'EDP est parabolique ; enfin, si la courbe est une hyperbole l'EDP est hyperbolique.

Exemple : en appliquant la méthode ci-dessus à l'équation(2.6), on obtient :

T² -- A.²X² = Cst (2.10)
qui est l'équation d'une courbe hyperbolique dans le plan(X, T).

1.4. EDPs du second d'ordre à plusieurs variables indépendantes:

Le principe de classification des EDPs reste le même quand on se trouve en présence de 3(ou 4) variable indépendantes. Par exemple, l'EDP (2.6) peut être généralisée à deux dimensions d'espace :

a²U a²u a²u

A.2 A² = 0 (2.11)

at2 x ax2 Y ay2

où A._x² et 4 sont les vitesses de propagation des ondes les directions x et y (elles sont égales dans un milieu isotrope). En utilisant la transformation exposée dans le paragraphe précédent. On obtient :

T² -- A._x²X²--A._Y²Y² = Cst (2.12)

1.5. Besoins en termes de conditions initiales et aux limites :

· EDPs hyperboliques :

Les EDPs hyperbolique abordées dans ce chapitre comprendrons en général une variable de temps et une (ou deux) d'espace. Elles seront de la forme (2.6). On cherche une solution de (2.6) en tout point d'un domaine de calcul SI = [x1,x2] et pour des dates t 0. Pour pouvoir déterminer la solution de (2.6) de façon unique, il faut connaître :

· la valeur de U en tout point du domaine de calcul à t = 0 (condition initiale) ;

· la valeur de U à tout date t 0 aux limites du domaine (condition aux limites) ;

Comme il y a des EDPs hyperbolique n'ont besoin que d'une condition à la limite.

· EDPs paraboliques :

La plupart des EDPs paraboliques de l'ingénieur sont d'ordre 1 par rapport au temps et d'ordre 2 par rapport à une (ou plusieurs) dimension(s) d'espace. C'est le cas de l'équation (2.7). Les conditions aux limites nécessaires à l'existence et l'unicité de la solution sont les mêmes pour que les EDPs hyperboliques.

· EDPs elliptique :

Les EDPs elliptiques étudiées dans ce qui suit impliquent 2 dimensions de l'espace et aucune de temps. Cela signifie que la solution ne dépend pas du temps. Dans ce cas, seules des conditions aux limites sont nécessaires.

2. Méthode des différences finies pour les EDP paraboliques :

Dans ce paragraphe nous avons aborderons la résolution numérique de problèmes d'équations aux dérivées partielles d'évolution par la méthode des différences finies. Nous étudierons un type de problèmes d'évolution du premier ordre en temps, dénommés également problèmes paraboliques.

Les équations intervenant dans ces problèmes sont constituées pour partie d'une combinaison de dérivées partielles par rapport à la variable temporelle dont nous détaillerons le traitement numérique et pour partie d'une combinaison de dérivées partielles par rapport à la variable spatiale ; le problème pouvant être posé dans un domaine , monodimensionnel, bidimensionnel ou tridimensionnel; pour simplifier l'exposé nous considérerons que le domaine est le segment [0, 1], le cas bi et tridimensionnel ne présentant pas de difficultés majeures.

2.1. Problèmes du premier ordre en temps : équation de la chaleur :

· Position du problème :

On considère une barre métallique de longueur unité; on suppose que cette barre est soumise à un apport de chaleur f (x, t ) par unité de longueur et de temps et

que, de plus, la température u (x, t) de la barre est maintenue à zéro à chacune de ses extrémités. On désigne par C la capacité thermique et par K le coefficient de diffusion de chaleur ; rappelons que le coefficient de diffusion de chaleur exprime qu'en tous les points de la barre les températures ont tendance à s'uniformiser. En supposant que la dimension transversale de la barre est négligeable par rapport à sa dimension longitudinale, la modélisation de ce problème conduit à déterminer u (x, t) en chaque point x E[0,1] et pour tout instant t E[0,T ], T représentant l'horizon d'intégration, solution de l'équation de la chaleur :

~ ~ ~~~~, ~~

~~ ~ ~ ~~~~~, ~~

~~~ ~ ~~~, ~~, ~ ~ ~0,1~ ; ~ ~ ~0, ~~

~ ~~0, ~~ ~ ~~1, ~~ ~ 0, ~ ~ ~0, ~~

~

u(x, 0) = uo(x), x E [0,1]

avec uo(x) température initiale de la barre.

Remarque :

Les conditions aux limites de type Dirichlet homogènes u (0, t) = u (1, t) = 0 correspondent, évidemment, au fait que les extrémités de la barre sont maintenues à une température nulle.

Remarque :

Dans le cas d'un matériel non homogène, le problème précédent est remplacé par :

~ ~~~~ ~~~~, ~~ ~~ ~~~~~ ~~~~, ~~

~~ ~ ~ ~~~, ~~, ~ ~ ~0,1~ ; ~ ~ ~0, ~~

~~ ~ ~

~

~ ~~0, ~~ ~ ~~1, ~~ ~ 0, ~ ~ ~0, ~~

~

~ ~ ~ ~0,1~
~~~, 0 ~ ~~~~~, Dans la suite, on supposera pour simplifier que C = 1,K = a (constante positive) ;

on ne considérera de plus que les conditions aux limites de Dirichlet homogènes, ce qui conduira à la formulation du problème modèle suivant :

~~~~~, ~~

~~ ~ ~ ~~~~~, ~~

~~~ ~ ~~~, ~~, ~ ~ ~0,1~, ~ ~ ~0, ~~

~ ~~0, ~~ ~ ~~1, ~~ ~ 0, ~ ~ ~0, ~~

~

u(x, 0) = uo(x), x E [0,1]

2.2.1. Schémas numériques de discrétisation par différences finies :

La discrétisation du problème précédent consiste à remplacer par une technique appropriée le problème continu par un système linéaire algébrique. L'approximation s'effectue par la manière suivante:

~

on divise l'intervalle [0,1] en (n + 1) intervalles de longueur h = .

n+1

on divise l'intervalle [0, T ] en M intervalles de temps k, tels que T = Mk. On pose xi = ih, pour i E f 1 , 2 , ..., n}, t_m = mk pour m E (0, 1, . . . , M} ; on écrit le problème modèle au point (xi, t_m ) avec 1 i n et m 0:

(au ~~ ~ ~ ~~~

~~~~ ~~~, ~~~ ~ f(xi,t_i)

~~~~~,~~~ ~~~~~~,~~~

On remplace et par des quotients différentiels faisant intervenir

~~ ~~~

les valeurs de u aux points xi, 0 i n + 1, et aux instants t_m, m > 0 . Exemple : on écrira :

a^2u(xi, tm)		~~~~~~, ~~~ ~ 2, ~~~ ~ ~~~~~~, ~~~ ~ + O(h^2) h^2	ax^e

	au(xi, tm)		u(xi, tm+i) -- u(xi, tin) = + O(k) k	et on obtient ainsi :
at

u(xi, tm+i) -- u(xi, tm)

~

~

~~ ~~~~~~, ~~~ ~ 2, ~~~ ~ ~~~~~~, ~~~

~ ~ O, ~~ ~ ~~~~, ~~~

~~

~ 1 ~ ~ ~ ~ ~~ ~ ~ 0

u(xo, tm) = u(xn+1, tm) = 0 u(xi, 0) = uo(xi)

~

~

Cette écriture permet de définir les schémas numériques définissant um r-r^, u(xi, t_m), en traitant de manière particulière le terme temporel. La façon la plus générale de définir les schémas aux différences finies est de considérer une combinaison convexe de l'équation de la chaleur aux instants t_m et t_m+1 ; plus

précisément, soit 0 E [0,1] un nombre réel ; écrivons l'approximation de (1 - 0) fois le problème modèle considéré à l'instant t_m et de 0 le même problème à l'instant t_m+1, ce qui permet de définir le 0 -- schema suivant :

uro_. -- 2ui^m + ur 1r _Er -- 24¹⁺¹ + uTV

ui um-- a ((1 -- 0) + 0u

k h² h²

= (1 -- 0)fi^m + 0fi^m+1 , 1 < i < n et m > 0

~ ug = u_n^p ₊₁ = 0, p = m et p = m + 1

~ ui^o = uo(ih)

avec f_i^p = f(xi, t_p), pour p = m ou p = m + 1

Pour 0 = 0, le schéma précédent définit un schéma explicite, c'est-à-dire que les valeurs um étant connues (par la condition initiale lorsque m = 0, et pour m > 0 par le calcul de l'itéré en temps précédent), les valeurs sont ur^l déterminées par la récurrence :

1

ur + ¹ = k f iⁿ' + ur + a (ur^. 1-- 24 ' + 4¹: 1) , 1 < i < n et m > 0 um = u^m = 0 o n+1
ui^o = uo(ih)

avec a =

ak h^2.

Si l'on définit les vecteurs de dimension n par :

F^m = (kf~ ^m, ... , kfnm)t U^m = (u^m, ,un^m)^t et la matrice B de dimension n par :

2 --1

--1 2

B =
·..

)

--1 2 --1

--1 2

le schéma explicite s'écrit matriciellement :

U^m+1 = F^m + (I -- aB)U^m

Pour 0 # 0, le schéma précédent définit un schéma implicite, c'est-à-dire qu'à chaque pas de temps, en supposant connu le vecteur U^m (toujours par la condition

ainsi que le schéma rétrograde implicite, également à trois niveaux, suivant :

~~ ~~~ ~ ~~ ~~~ ~~~ ~ ~~ ~~~ ~ ~~~~

~

~~~~

~ ~ ~~

~ ~ ~

~ 2 ~~

1

~ ~~ ~ ~ ~~~~

~~~ ~ 0

~ ~ ~~ ~ ~ ~~~~~~

~ ~~ ~

~ ~ ~ ~ ~~ ~ ~ 0

initiale lorsque m = 0 et pour m > 0 par le calcul de l'itéré en temps précédent), le vecteur U^m+1 est obtenu en résolvant le système linéaire suivant:

(I -- a0B)U^m+1 = (1 -- 6)F^m + OF' + (I -- a(1 -- 6)B)U^m

Notons immédiatement que, pour 0 E0, 1], la matrice (I -- a0B) est inversible ; dans le cas du problème monodimensionnel en espace, la résolution de ce système s'effectue par une simple adaptation de la méthode de Gauss, (pour l'adaptation de cet algorithme au cas des matrices tridiagonales). Pour des valeurs particulières du paramètre 0, on retrouve des schémas numériques particuliers ; ainsi pour 0 = 1, on obtient le schéma purement implicite qui s'écrit matriciellement :

(I -- aB)U^m+1 = F^m+1 + U^m

alors que pour 0 = 0,5, on obtient le schéma implicite de Crank-Nicolson:

1 a

(I --^a₂ 13) U^m+1 = 2 (F^m + F^m+1) + -- 2 /3) U^m

Remarque :

Il existe d'autres types de schémas numériques pour résoudre le problème modèle, par exemple le schéma saute-mouton défini par :

~ ~

~~ ~~~ ~ ~~

~ ~ ~~ ~

2 ~ ~ ~~~~

~ ~ 2 ~ ~ ~~~~

1 ~ ~ ~ ~ ~~ ~ ~ 0

~~

~ ~ ~~~~

~ ~ 0

~ ~~

~

u° = uo(ih)

C'est un schéma à trois niveaux explicite qui, bien que plus précis que les schémas précédents, n'est pas intéressant sur le plan numérique car il est instable, c'est-à-dire sensible à la propagation des erreurs systématiques d'approximation, de chute ou de troncature. Citons également le schéma explicite à trois niveaux proposé par Du Fort et Frankel:

~~~

~ 1 ~

~3 ~~~~ ~ ~ ~~~~

~~~ ~ 2

2 ~~ ~ ~ ~~~~

2 ~~ ~~~ ~ 2 ~ ~ 1 ~ ~~

~ ~ ~~

1 ~ ~ ~ ~ ~~ ~ ~ 0

~ ~~ ~~~ ~ ~~~~

= 1.10(ih)

~~~ ~ 0

~ ~ ~~ ~

2.2.2. Erreur de troncature, consistance et ordre d'un schéma :

~

Pour chacun des schémas précédents, on remplace dans le schéma le terme ~~ par

u(xi, t_m) et l'on définit l'erreur de troncature, notée, comme la différence entre le premier et le second membre de ces quantités. Ainsi pour le 0 -- schéma considéré au paragraphe précédent, on a l'expression suivante de l'erreur de troncature :

~~~~, ~~~~~ ~ ~~~~, ~~~ ~ ~ ~~1 ~ ~~ ~~~~~~, ~~~ ~ 2, ~~~ ~ ~~~~~~, ~~~

~~ ~ ~ ~

~~

~~~~~~, ~~~~~ ~ 2, ~~~~~ ~ ~~~~~~, ~~~~~

~ ~ ~ ~ ~1 ~ ~~~~~~, ~~~ ~ ~~~~~, ~~~

~~

Remarque :

On définit l'erreur de troncature associée au schéma saute-mouton, au schéma rétrograde et au schéma de Du Fort et Frankel de manière analogue.

On adapte à la situation des problèmes d'évolution les notions d'erreur de troncature, de consistance et d'ordre de schéma comme suit :

Définition 1 : On appelle erreur de troncature E, la quantité définie par : E = Max(E1 ⁿ, 1 < i < n,1 < m < M}

Définition 2 :

Un schéma numérique est consistant si E --> 0 lorsque h --> 0 et k --> 0. Définition 3 :

Un schéma est précis à l'ordre p en espace et à l'ordre q en temps, s'il existe une constante C', indépendante de h et k telle que |E| Cⁱ(h^P + k ^q)

On a le résultat suivant:

Théorème 1 :

On suppose que la solution du problème modèle est suffisamment régulière ; alors pour 0 E [0, 1], 0 # 0,5, le 0-schéma est d'ordre 2 en espace et d'ordre 1 en

temps. Pour 0 = 0,5, le schéma de Crank-Nicolson est d'ordre 2 en espace et en temps.

Remarque : Il résulte du résultat précédent que le schéma explicite et le schéma purement implicite sont d'ordre 2 en espace et d'ordre 1 en temps.

Remarque :

On vérifie que le schéma saute-mouton et le schéma rétrograde sont d'ordre 2 en espace et en temps. L'étude du schéma de Du Fort et Frankel est plus délicate; en

effet, ce dernier n'est consistant que si le rapport kh --> 0 lorsque h --> 0 et k --> 0.
De plus, si l'on veut que le schéma de Du Fort et Frankel soit d'ordre 2 en espace, il faut prendre k = 0(h²).

2.2.3. Stabilité des schémas numériques :

La stabilité d'un schéma numérique pour l'approximation de la solution d'une équation aux dérivées partielles est une question importante et, apparemment, nouvelle ; lorsque le pas de discrétisation tendait vers zéro, le système linéaire pouvait devenir mal conditionné, c'est-à-dire que la solution calculée du système discret pouvait être dénaturée ; en effet, la question qui se pose peut être formulée de la façon suivante : compte tenu d'une part des erreurs systématiques introduites lors de la construction du schéma d'approximation et d'autre part compte tenu également de la mauvaise représentation des nombres en machine qui induisent des erreurs d'arrondi ou de troncature, est-ce que l'erreur va rester bornée et tendre vers zéro lorsque les paramètres de discrétisation h et k vont tendre vers zéro ?

Il existe de nombreuses méthodes pour étudier la stabilité d'un schéma numérique ; nous en étudierons deux, connues sous le nom de :

méthode matricielle ;

méthode de Fourier ou méthode du coefficient d'amplification.

La méthode matricielle : consiste à considérer qu'un schéma est stable si, lorsqu'une erreur est introduite sur les conditions initiales (par exemple, erreur d'arrondi ou de troncature), l'erreur sur la solution calculée à un instant fixé t n'est pas amplifiée ; l'idéal est, bien entendu, qu'elle ne soit pas amplifiée du tout et que

supi,m(Iui^m - u(xi, t_m)1) reste borné quels que soient h et k. L'étude de la stabilité d'un schéma revient donc en fait à celle de l'amplification d'une perturbation sur la condition initiale U^°; on va vérifier que l'erreur sur la solution cherchée satisfait le schéma numérique homogène associé, ce qui va permettre de dégager des critères de stabilité. A ce stade, on peut donc donner la définition intuitive suivante de la stabilité :

Definition :

Un schéma numérique est stable si la solution du schéma homogène associé est bornée en tout point xi et t_m quels que soient h et k.

Pour illustrer cette approche intuitive, considérons le 0 - schéma dont la formulation est la suivante :

(I - a0B)U^m+1 = (1 - 0)F^m + 0F^m+1 + (I - a(1 - 0)B)U^m

et soit U^m+1 la solution exacte du schéma numérique et V^m+~ celle du schéma perturbé : on peut traduire l'imprécision par la relation suivante :

V^m = U^m + E^m, Vm > 0

avec E^m erreur au pas m.

En remplaçant U^m+1 par V^m+~ dans le schéma précédent, on obtient :

(I - a0B)E^m+1 = (I - a(1 - 0)B)E^m

et on vérifie ainsi que le terme d'erreur est solution du schéma numérique homogène associé. Cette relation peut encore s'écrire :

E^m+1 = (I - a0B)^-1(I - a(1 - 0)B)E^m = PE^m =
· = Pm+1E° avec P = (I - a0B)^-1(I - a(1 - 0)B).

Le schéma numérique sera donc stable si la norme de la matrice P^mreste bornée lorsque m -> 00. La matrice B étant symétrique, on vérifie aisément que la matrice P est normale (c'est-à-dire que P^tP = PP^t ) ; on peut alors exprimer la relation précédente dans la base des vecteurs propres de la matrice P ; soit A la matrice des valeurs propres de la matrice P. La relation E^m+1 = P^m+1E^°, s'écrit alors E~^m+~ = A^m+1E^°, où E~^m+~ est l'écriture de E^m+1 dans la base des vecteurs propres, et une condition suffisante de stabilité est donc :

|uk(P) | < 1,V k E (1,...,n )

avec uk(P) valeurs propres de la matrice P, c'est-à-dire les coefficients diagonaux de la matrice .4.

On vérifie facilement que :

1 -- a(1 -- 0)A.k(B)

uk(P) 1 + a0A.k(B) '

= V k E (1,...,n )

avec A.k(B) valeurs propres de la matrice B; on vérifie alors que la

condition - 1 < uk(P) < 1, V k E (1, ... , n ), se traduit par la condition suivante :

ak a = < h^2

1 , 0 E [0,1] 2(1 -- 20)

qui représente donc la condition de stabilité. On peut résumer l'étude précédente par le théorème suivant :

Théorème 2 :

Le 0-schéma est inconditionnellement stable si 0 > 0,5 ; si 0 < 0,5, ce schéma est stable si :

ak

a = <

h²

1

2(1 -- 20)

Corollaire 1 :

Le schéma purement implicite (0 = 1) et le schéma de Crank-Nicolson (0 = 0,5) sont inconditionnellement stables. Le schéma purement explicite (0 = 0) est stable à condition que :

1

2

ak

a = <

h²

Remarque :

Le schéma saute-mouton est instable ; le schéma de Du Fort et Frankel et le schéma rétrograde sont inconditionnellement stables.

La méthode matricielle présentée ci-dessus est applicable dans le cas où la matrice B est symétrique. Or il existe des situations où cette matrice ne vérifie pas cette condition ; c'est le cas par exemple du problème de convection-diffusion d'évolution. Dans ce cas, le calcul du coefficient d'amplification permet de s'affranchir de cette

contrainte et fournit une méthode d'analyse plus générale de la stabilité. La définition de la stabilité étant inchangée, considérons toujours le 0-schéma homogène :

ur^l -- -- a[(1 -- ~~~~~~~

~ -- 2ur + ~~~~

~ ~

~~~~~~~

~~~ -- 24¹⁺¹ + ~~~~

~~~~~ = 0,1 < i < n et m > 0

que l'on peut encore formellement écrire en repassant à la variable d'espace continue x et en mettant dans le premier membre les termes à l'instant discret (m + 1) et dans le second membre ceux à l'instant discret m :

u^m+1(x) -- a0(u^m+1(x + h) -- 2u^m+1(x) + u^m+1(x -- h))

= u^m(x) + al1 -- 0) ~~~~~

~~~~~ + h) -- 24¹⁺¹(x) + ~~~~

~~~~~ -- h)) , m > 0

Considérons la transformée de Fourier du schéma ainsi écrit ; on rappelle que cette transformation est définie par :

~~

1

~~~~~ ~ _v2~ ~ ~~~~~~~~~~~~~ ~~, ~~ ~ ~1

~~

On vérifie aisément par un calcul direct que :

~~

1

v2~ ~ ~~~ ~ ~~~~~~~~~~~ ~~ ~ ~~~~~~~~~~~~~~ _00

et le schéma discret en temps s'écrit alors :

It^m+1(0 = c(a)u^m(

où c() est le coefficient d'amplification défini par :

1 + 2a(1 -- 0)(cos(~h) -- 1)

~~~~ ~1 -- 2a0(cos(~h) -- 1)

~

1 -- 4a(1 -- 0)sin2 ~~~ 2 ~

1 + 4a0sin² (^Q1₂ )

=1

asin2 ~~~

2 ~

1 + a0sin2 ~~~

2 ~

Comme précédemment, une condition suffisante de stabilité est donnée par --1 c(0 1 ; or on a évidemment c() < 1; la vérification de la seconde
inégalité conduit à la condition suivante :

ak a= < h^2

1 , ~ ~ ~0,1~ 21 ~ 2

qui correspond à la condition obtenue par la méthode matricielle ; on retrouve donc les résultats obtenus et énoncés dans le théorème 2 et le corollaire 1.

2.2.4. Convergence des schémas :

Soit e^m , le vecteur erreur au temps discret m et défini par :

e^m = (er, ... , ei ^m, ... , ein^t

avec :

~~ ~ ~~ ~ -- u(xi,t_m), 1 i n,m 0

Définition :

Un schéma numérique est convergent si :

|~T -- u(xi, t_m)| --) 0 lorsque h --) 0 et k --) 0 ,Vi e (1,...,n),Vm 0.

Considérons le 0- schéma ; en considérant d'une part l'écriture du schéma et d'autre part celle de l'erreur de troncature et en soustrayant membre à membre, on vérifie aisément que les composantes de l'erreur satisfont le schéma numérique, avec un second membre égal à --Elⁿ, soit :

~~~ ~ ~~ ~

~ ~~ ~ ~ 2 ~ ~ ~~~~

~

~ ~~1 ~ ~~ ~~~~

~ ~~

~ ~~~

~~~~

~~~ ~ 2 ~~~ ~ ~~~~

~~ ~ ~~~ ~, 1 ~ ~ ~ ~ ~~ ~ ~ 0

~ ~~

~ ~ ~ ~ ~ 1

~ ~~ ~ ~~

~ ~~~~ ~ 0, ~ ~ ~

~

soit matriciellement :

(I -- a0B)e^m+1 = E^m + (I -- a(1 -- 0)B)e^m

où les composantes du vecteur E^msont égales à l'erreur de troncature, au signe près ; en multipliant scalairement par e^m+1, il vient :

((I -- a0B)e^m+1, e^m+1) = (E^m, e^m+1) + ((I -- a(1 -- 0)B)e^m, eⁿ⁺¹),Vm 0 Or, la matrice(I -- a0B) étant symétrique définie positive, on peut minorer le second membre; si, de plus, on note par p , la norme matricielle de (I -- a(1 -- 6)B) , on obtient finalement :

eP ~ = 0

Ile^m+1112 ^filEm12 _+pllem112,vm 0

avec /3 inverse de la plus petite valeur propre de la matrice (I -- a0B).

En posantT^m = Ile^m112,Vm 0, on obtient finalement :

T^m+1 /311E^m112 + pT^m, Vm 0

Supposons, de plus, que le schéma vérifie la condition de stabilité, ce qui a pour conséquence que la quantité T^m = Ile^m112, Vm 0, est borné. Pour pouvoir obtenir une majoration de la norme de l'erreur, on utilise le lemme suivant dont le résultat s'obtient par récurrence :

Lemme 1 : Si l'on a une suite qui vérifie la relation de récurrence :

T^m+1 pT^m + a,Vⁿ 0,Vm 0, p # 1

alors :

p^m -- 1

T^m < pT^° + (5 IlE^m112,Vm 0

p -- 1

En appliquant ce résultat à notre situation, compte tenu du fait que E^° = 0 , on obtient :

~~~ ~ ~~~~~ ~ ~ ~~ ~ 1

~ ~ 1 , ~~ ~ 0

Si, de plus, le schémaa numériquee est consistant, l'inégalitée précédentee nous permet de vérifierr qu'ill est convergent. On peut donc énoncerr le résultatt suivant ^:Théorèmee 3 : Si le 0 - schémaa est stable et consistant alors ilt est convergent.

Le résultatt précédentt exprime que la stabilitée et la consistance constituent une condition suffisante pour la convergence. En fait le théorèmee de Lax( le théorèmee de Lax prévoitt que dans un problèmee bien posé,, et avec un schémaa numériquee consistant, la stabilitée est une condition nécessairee et suffisante pour la convergence), nous donne le résultatt suivant :

Théorèmee 4 :

La stabilitée est une condition nécessairee et suffisante pour qu'unn schémaa numériquee soit convergent vers la solution du problème,, si ce dernier est discrétisée par un schémaa consistant.