85
  MODELISATION Introduction aux équations aux dérivées partielles et à leurs résolutions numériques FN CRES

cours de modélisation

Embed Size (px)

Citation preview

Page 1: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 1/85

 

MODELISATION

Introduction aux équations aux dérivéespartielles et à leurs résolutions numériques

FN CRES

Page 2: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 2/85

Modélisation

FN CRES

2

 

SOMMAIRE

AVANT-PROPOS ......................................................................................................................... 4 

1.  INTRODUCTION ................................................................................................................... 5 

2.  EXEMPLE D'ETABLISSEMENT D’UNE EQUATION MECANISTE ............................................6 

2.1. Equation de diffusion 1D ............................................................................................. 6

2.2. Interprétation d'une équation aux dérivées partielles ................................................... 8

2.3. Equation de la Diffusion 3D ......................................................................................10

2.4. Conditions limites et condition initiale ......................................................................11

2.5. Régime permanent ..................................................................................................... 12

2.6. Compléments à l’équation de base............................................................................. 12

2.6.1.   Injection (ou soutirage) ....................................................................................12 

2.6.2.  Cinétique propre du produit............................................................................. 12 

2.6.3.  Convection........................................................................................................ 13 

2.6.4.  Cas où K est variable (dans toutes les directions et en tous points)................ 13 

2.6.5.   Autres problématiques...................................................................................... 15  

3.  CLASSIFICATION DES EQUATIONS AUX DERIVEES PARTIELLES D’ORDRE 2 ...................16 3.1. Classification par les coniques................................................................................... 16

3.1.1.  Classification par discriminant et par valeurs propres ................................... 16  

3.1.2.   Application aux EDP d’ordre 2 ....................................................................... 16  

3.2. Equation stationnaire et équation d’évolution ........................................................... 17

4.  RESOLUTION DES EDP - PRINCIPE DES METHODES NUMERIQUES .................................. 18 

5.  LA METHODE DES DIFFERENCES FINIES ........................................................................... 19 5.1. Principe de la méthode............................................................................................... 19

5.2. Cas d’une EDP elliptique (stationnaire)..................................................................... 20

5.2.1.  Cas de conditions de Dirichlet ......................................................................... 20 

5.2.2.  Cas de conditions de Neumann homogènes..................................................... 26  

5.2.3.  Cas de conditions de Neumann non homogènes .............................................. 29 

5.3. Notions sur les erreurs de la méthode ........................................................................ 31

5.3.1.  Présentation ..................................................................................................... 31 

5.3.2.   La consistance.................................................................................................. 32 

5.3.3.   La stabilité........................................................................................................33  

5.3.4.   La convergence ................................................................................................ 33 

5.4. Cas d’une EDP parabolique (évolutive).....................................................................335.4.1.  Position du problème et méthode.....................................................................33 

5.4.2.   Les schémas types............................................................................................. 34  

5.4.3.   Le schéma explicite .......................................................................................... 36  

5.4.4.   Le schéma implicite .......................................................................................... 36  

5.4.5.   Le schéma explicite à 2 pas..............................................................................38 

5.4.6.   Le schéma pondéré........................................................................................... 38 

5.4.7.  Cas d'un schéma explicite stable...................................................................... 40 

5.4.8.  Cas d'un problème parabolique 2D .................................................................40 

5.5. Cas d'une EDP hyperbolique (évolutive) ...................................................................41

5.5.1.  Schémas explicites............................................................................................ 42 

5.5.2.  Schémas implicites ........................................................................................... 43 5.6. Maillages particuliers ................................................................................................. 45

Page 3: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 3/85

Modélisation

FN CRES

3

5.6.1.   Maillage irrégulier en x ...................................................................................45 

5.6.2.   Maillage différent en x et en y.......................................................................... 46  

5.7. Traitement des termes non linéaires........................................................................... 47

6.  LA METHODE DES ELEMENTS FINIS ..................................................................................48 

6.1. Comparaison Différences finies – Eléments finis...................................................... 48

6.2. Equivalence de problèmes .........................................................................................486.2.1.  Objet ................................................................................................................. 48 

6.2.2.   Equivalence entre système matriciel et problème de minimisation .................48 

6.2.3.   Equivalence avec un problème de minimum.................................................... 49 

6.2.4.   Récapitulatif ..................................................................................................... 50 

6.2.5.   Application aux EDP........................................................................................ 50 

6.2.6.   Approximation interne...................................................................................... 51  

6.3. Construction pratique du problème variationnel........................................................ 52

6.3.1.  Cas d’une équation différentielle ..................................................................... 52 

6.3.2.  Cas d’une EDP................................................................................................. 54 

6.4. Mise en œuvre de la méthode des éléments finis ....................................................... 55

6.5. Cas d’une équation différentielle - Fonction linéaire par morceaux.......................... 556.5.1.  Cas de condition de Dirichlet homogène ......................................................... 55 

6.5.2.  Cas de condition de Neumann homogène ........................................................ 60 

6.5.3.  Cas de condition de Neumann non homogène ................................................. 61 

6.5.4.  Cas de condition de Dirichlet non homogène .................................................. 62 

6.6. Cas d’une équation différentielle - Fonction parabolique par morceaux.................. 63

6.7. Cas d’une EDP – Approximation linéaire par élément.............................................. 70

6.8. Cas d’une EDP – Approximation quadratique par élément ....................................... 77

6.9. Extension de la méthode ............................................................................................78

6.10. Extension aux problèmes évolutifs ............................................................................ 78

7.  RESOLUTION DES SYSTEMES D’EQUATIONS LINEAIRES................................................... 79 

7.1. Introduction................................................................................................................79

7.2. Méthodes directes ...................................................................................................... 79

7.2.1.   Méthode de Gauss-Jordan ou méthode du pivot .............................................. 80 

7.2.2.   Méthode de Gauss (ou triangularisation) ........................................................ 81 

7.3. Méthode du double balayage pour les matrices tridiagonales (Cholesky)................. 82

7.4. Méthodes itératives .................................................................................................... 83

7.4.1.  Principe ............................................................................................................83 

7.4.2.   Méthode de Jacobi ........................................................................................... 84 

7.4.3.   Méthode de Gauss -Seidel................................................................................ 84 

7.4.4.  Facteur de relaxation.......................................................................................85 

Page 4: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 4/85

Modélisation

FN CRES

4

AVANT-PROPOS

Ce polycopié est une introduction aux équations aux dérivées partielles et à leur résolution,tout ceci ayant comme objectif la modélisation mathématique du monde qui nous entoure, ou,

 pour rester plus modeste, la modélisation des problèmes courants rencontrés par l'ingénieur.

Le chapitre 1 situe succinctement la modélisation mathématique parmi les modélisations les

 plus courantes.

Ces équations ne sont pas toujours bien abordées par les étudiants qui y voient une écriture

ésotérique et parfois incompréhensible. Le chapitre 2 tente de démystifier ces équations en

 présentant un exemple simple et baser sur "le bon sens" dans lequel on aboutit à l'équation de

la diffusion. L'objectif est alors de ne pas perdre le lecteur avec une approche trop

mathématique susceptible de le décourager. Cependant, il est important que le lecteur se

familiarise rapidement avec les symboles mathématiques utilisés.

Ce n'est pas la manière "officielle" de présenter cette équation, mais cette démonstration

 présente l'immense avantage de pouvoir être suivie avec un niveau de terminale scientifique

(du moins, je l'espère !). On arrive petit à petit à l'écriture la plus complète de l'équation de la

diffusion.

Le chapitre 3 n'est pas le plus important, mais il introduit la classification des équations aux

dérivées partielles, élément structurant pour la suite.

A partir du chapitre 4, on aborde la résolution des équations aux dérivées partielles par 2

méthodes numériques :

•  La méthode des différences finies, qui est une méthode simple, et que lecteur pourra

suivre aisément

•  La méthode des éléments finis, qui est d'un abord beaucoup plus complexe. Le lecteur 

s'attachera à comprendre le principe et à traiter les exemples complets qui y sont

développés.

Enfin, le dernier chapitre sur la résolution des systèmes d'équations linéaires, aboutissementdes 2 méthodes précédentes, est là à titre informatif. En effet, la mise en œuvre de ces

méthodes nécessite généralement des connaissances informatiques, notamment de

 programmation, ce qui déborde du sujet abordé dans ce document.

Page 5: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 5/85

Modélisation

FN CRES

5

 

1. INTRODUCTION

La modélisation d’un phénomène est une démarche visant à représenter par un moyen adéquat

le comportement de ce phénomène. Dans les sciences de l'ingénieur, la modélisation permet

de comprendre les variables qui influencent ce comportement, afin de dimensionner desouvrages, d'anticiper son évolution, de simuler des situations à venir.

La modélisation peut être abordée de différentes façons. On peut proposer la classification

sommaire suivante :

•  Les modèles réduits

Qui ne connaît pas les souffleries où sont testés les modèles réduits d’avion ? Le modèle

réduit permet de rendre compte du comportement d’un objet soumis à différentes

contraintes sans avoir à construire cet objet dans sa taille normale. La théorie des

similitudes permet alors, à partir du comportement du modèle réduit, de conclure sur le

comportement de l’objet réel.

•  Les modèles analogiques

Ils permettent de représenter un phénomène à partir d’une analogie avec un autre plus

facile à élaborer. Par exemple, le comportement d’une nappe d’eau dans le sol peut être

abordé par une analogie avec le potentiel électrique d'une plaque métallique.

•  Les modèles mathématiques

Ces modèles sont les plus courants actuellement, suite à la montée en puissance des

ordinateurs et de leur capacité à calculer vite. Ils sont basés sur la mise en équation

mathématique du phénomène à étudier. Ce sont ces modèles qui vont nous intéresser pour ce cours. Là aussi, on peut tenter une classification sommaire.

-  Les modèles empiriques

Il s’agit d’identifier les variables qui interviennent à priori dans un phénomène

  physique et de les relier par une équation à partir d’une série d’observations. Cette

équation n’a parfois rien de physique, mais représente bien le « nuage de points ». Elle

est totalement dépendante de l’échantillon qui a servi au calage.

-  Les modèles conceptuels

Ils abordent la représentation d’un phénomène complexe à partir d’un autre beaucoup

  plus simple à étudier. Par exemple, en hydrologie, on conçoit souvent lefonctionnement d’un bassin versant (en termes de production d’un débit d’eau) comme

celui d’un réservoir, objet dont le remplissage et/ou la vidange se mettent facilement

sous forme d’équations.

-  Les modèles mécanistes

La mécanique, en tant que science, est à la base de la représentation du phénomène.

On aboutit généralement à un type d’équations dites aux dérivées partielles, qu’il

s’agit ensuite de résoudre.

C’est ce dernier type de modélisation auquel ce cours se consacre. On verra successivement

comment on aboutit à des équations aux dérivées partielles à travers un exemple, et deux

méthodes classiques pour les résoudre.

Page 6: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 6/85

Modélisation

FN CRES

6

 

2. EXEMPLE D'ETABLISSEMENT D’UNE EQUATION MECANISTEOn va partir d’un exemple suffisamment simple pour être compréhensible quelque soit

l’origine scientifique du lecteur : le comportement d’un produit – par exemple un polluant – 

dans de l’eau.

Il s’agit d’un problème dit de diffusion.

2.1. Equation de diffusion 1DOn considère un parallélépipède, de section S, constitué de matière homogène immobile (de

l'eau par exemple – on verra plus loin le cas de l'eau en mouvement), ayant une concentration

C1 d'un produit sur sa face gauche et C2 sur sa face droite.

On peut faire comme hypothèse que la quantité

de matière M –issue du produit en question– qui

franchit une section S du parallélépipède, c'est-

à-dire qui circule par unité de longueur sur l'axe

des x, de l'avant vers l'arrière pendant le temps

Δt est:

-   proportionnel à la section S

-   proportionnel à la différence C1-C2 

-   proportionnel à Δt

-  inversement proportionnel à Δx. en effet,

 plus Δx est petit et plus la quantité de

matière devant franchir la section S

 pendant l'intervalle de temps Δt sera importante pour passer de la concentration C1 et

C2.

soit tx

CCSK M 21

x ΔΔ

−=  

 Remarque : cette relation peut se vérifier expérimentalement.

 La   forme de cette relation s'applique à différents domaines scientifiques en fonction de la

variable C et porte différents noms selon le domaine concerné (loi de Fick en Génie des

Procédés, loi de Darcy en hydrogéologie, loi de Fourier en thermique …)

On admettra comme convention que M est positif si C1 > C2 et que le coefficient K x est

constant le long de l’axe Ox (on étudiera le cas de ce coefficient variable plus loin)

Cette équation très simple, admise comme hypothèse, est à la base de l'équation de la

diffusion. Le lecteur est invité à bien la mémoriser et toujours se rappeler ce point de départ.

On considère maintenant un volume de matière homogène découpé en parallélépipèdes de

longueur  Δx. On considère des parallélépipèdes suffisamment petits pour que la

concentration, au sein de chaque parallélépipède, y soit considérée comme constante.

Δx

C1 

C2 

S

Page 7: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 7/85

Modélisation

FN CRES

7

Faisons un bilan de produit au niveau de la tranche i selon l'axe des x : on regarde ce qui

rentre et ce qui sort de cette tranche pendant l’intervalle de temps Δt afin de connaître la

concentration dans cette tranche. On suppose donc qu'il n'y a pas d'échange de matières dans

les directions Oy et Oz.

Ci-1 Ci Ci+1

1 2

x

z

y

Δx

Tranches

i-1 i i+1 

La quantité de matière qui franchit la face 1, comptée positivement dans le sens de l'axe Ox,

s'écrit :

tx

CCSK M i1i

x1i ΔΔ

−= −  

La quantité de matière qui franchit la face 2, comptée positivement dans le sens de l'axe Ox,

s'écrit :

t

x

CCSK M 1ii

x2i Δ

Δ

−= +  

Si on considère que le coefficient K x est constant le long de l'axe Ox, le bilan dans la tranche i

s'écrit :

accumulation (Mi) = ce qui entre (Mi1) – ce qui sort (Mi2)

tx

CC2CSK MMM 1ii1i

x2i1ii ΔΔ

+−=−= +−  

Ecrivons la variation de concentration dans la tranche i, en divisant par son volume S×Δx :

t

x

CC2CK t

xxS

CC2CSK 

xS

MC

2

1ii1ix

1ii1ix

i Δ

Δ

+−=Δ

Δ×Δ×

+−=

Δ×

=Δ +−+−  

Soit2

1ii1ix

x

CC2CK 

t

C

Δ+−

=ΔΔ +− (équation A)

Regardons de plus près les deux termes de cette équation :

  Quand Δt → 0, on peut écrire le premier terme de cette équation sous forme

différentiellet

C

∂∂

 

  A quoi correspond le deuxième terme de cette équation ?

Page 8: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 8/85

Modélisation

FN CRES

8

Considérons la courbe ci-dessous et sa dérivéex

C

∂∂

au point d'abscisse x. Cette dérivée

correspond à la tangente à la courbe au point d'abscisse x.

On peut approximer cette tangente

  par la pente de la corde entre les points (xi ,Ci) et (xi+1 ,Ci+1).

Cette pente est égale àx

CC i1i

Δ−+ ;

 plus Δx est petit, et plus

l'approximation est valide.

Autrement dit, le termex

CC i1i

Δ

−+ est

une approximation de la dérivée

 premièreix

C⎟

 ⎠

 ⎞⎜⎝ 

⎛ ∂∂

à l'abscisse xi quand Δx→ 0.

De la même façon, on peut approximer la dérivée premièreix

C⎟

 ⎠

 ⎞⎜⎝ 

⎛ ∂∂

à l'abscisse xi par la pente

de la corde entre les points (xi ,Ci) et (xi-1 ,Ci-1) :x

CC

x

C 1ii

i Δ−

≈⎟ ⎠

 ⎞⎜⎝ 

⎛ ∂∂ − quand Δx→ 0.

On peut maintenant avoir une approximation de la dérivée seconde en x i , en dérivant 2 fois

C(x) par rapport à x, et en appliquant les 2 façons vues au-dessus:

2

1ii1i

2

1ii

2

i1i

i1i

i1i

ii

2

2

x

CC2C

x

CC

x

CC

x

C

xx

C

xx

CC

xx

C

xx

C

Δ+−

=Δ−

−Δ

−=⎟

 ⎠

 ⎞⎜⎝ 

⎛ Δ∂

∂−⎟

 ⎠

 ⎞⎜⎝ 

⎛ Δ∂

∂=⎟

 ⎠

 ⎞⎜⎝ 

⎛ Δ

−∂∂

≈⎟ ⎠

 ⎞⎜⎝ 

⎛ ∂∂

∂∂

=⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂ +−−+

+

+

 

Autrement dit, le deuxième terme de l'équation A correspond à l'approximation de la dérivée

seconde de la fonction C; quand Δx → 0, nous avons

i

2

2

2

1ii1i

x

C

x

CC2C⎟⎟

 ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂

≈Δ

+− +−  

L'équation A est donc une expression de

2

2

x x

CK 

t

C

∂=

∂ 

et tend vers cette dernière quand Δx → 0 et Δt → 0

L'équation encadrée est dite de la diffusion 1D (1D pour une dimension d'espace).

C’est une équation aux dérivées partielles (EDP) car C est fonction d’au moins 2 variables (x

et t) et elle fait intervenir les dérivées de C par rapport à ces 2 variables.

2.2. Interprétation d'une équation aux dérivées partielles

Comment "lire" l'équation précédente ?

Ci

x

Xi Xi+1Xi-1

Δx Δx

Ci-1

Ci+1

C

tangente

Approximation

de la tangente

Page 9: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 9/85

Modélisation

FN CRES

9

La dérivée premièret

C

∂∂

correspond à la variation de concentration dans le temps, c'est-à-dire

de quelle quantité C varie dans le temps.

La dérivée seconde2

2

x

C

∂correspond à la variation de la variation (dérivée de la dérivée

 première) de concentration selon l'axe des x, autrement dit, ce terme s'intéresse à la façon dont

la dérivée première de la concentration varie en x.

L'équation de la diffusion dit que la variation de concentration au cours du temps est

 proportionnelle à la variation de la variation de concentration en espace.

Interprétons ce qui vient d'être dit :

x

Ci-ε 

C i+ε 

Xi+ε XiX i-εXi+εXi-ε  Xi

Cas 1 Cas 2

C

Sur le graphique ci-dessus, on considère un petit déplacement à gauche et à droite de X etl'allure de la solution entre Xi-ε et Xi+ε. Nous avons une valeur (approximative) des dérivées

 premières à droite et à gauche de Xi à travers les pentes des cordes entre Xi-ε et Xi d'une part,

et Xi et Xi+ε d'autre part.

Il est important de remarquer que la variation de concentration entre X i-ε et Xi+ε est ma même

dans les 2 cas.

Dans le cas 1, les pentes sont voisines. L'équation dit donc que dans l'intervalle de temps Δt,

la concentration en xi ne changera guère (donct

C

∂∂

sera faible) puisque le terme2

2

x

C

∂∂

est voisin

de 0.Dans le cas 2, les pentes sont très différentes. L'équation de la diffusion dit alors que la

concentration en xi va beaucoup changer (t

C

∂∂

sera important) puisque le terme2

2

x

C

∂∂

est loin

d'être négligeable).

Ainsi, l'équation de la diffusion précise que la variation de concentration dans le temps ne

dépend pas de la différence de concentration entre 2 points, mais dépend de la différence des

 pentes entre ces 2 points (donc de l'allure de la solution et non de la valeur de la solution).

L'interprétation physique est simple : les valeurs de pentes représentent les flux de part etd'autres du point xi (soit des quantités de polluant qui circulent à gauche et à droite de xi , et

Page 10: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 10/85

Modélisation

FN CRES

10

c'est bien ce qu'on a écrit au tout début : on a fait un bilan de ce qui entrait et sortait de la

tranche i).

-  Si les flux de part et d'autre de x i sont voisins, le bilan est proche de zéro et la

concentration en xi ne varie guère dans le temps.

-  Si les flux sont différents, le bilan ne sera pas équilibré au niveau de x i et la

concentration changera dans le temps

2.3. Equation de la Diffusion 3DEn supposant qu'il y a également des échanges dans les directions Oy et Oz.

On peut réaliser ce même bilan de matière selon les deux autres axes y et z de l'espace, pour 

avoir un bilan total au niveau de la tranche C i,j,k .(indice i pour l’axe Ox, j pour l’axe Oy et k 

 pour l’axe Oz)

Supposons dans un premier temps que le coefficient de proportionnalité est constant dans

toutes les directions de l'espace : K x = K y = K z = K 

L'accumulation Mi,j,k de matière dans la tranche (i,j,k) s'écrit :

tz

CC2CSK 

ty

CC2CSK 

tx

CC2CSK 

)MM()MM()MM(M

1k , j,ik , j,i1k , j,i

zz

k ,1 j,ik , j,ik ,1 j,i

yy

k , j,1ik , j,ik , j,1i

xx

2k , j,2i1k , j,ik ,2 j,ik ,1 j,ik , j,2ik , j,1ik , j,i

ΔΔ

+−+

ΔΔ

+−+

ΔΔ

+−=

−+−+−=

+−

+−

+−

 

Supposons dans un premier temps que le coefficient de proportionnalité est constant danstoutes les directions de l'espace : K x = K y = K z = K 

⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

Δ

+−+

Δ

+−+

Δ

+−Δ= +−+−+−

z

CC2CS

y

CC2CS

x

CC2CStK M

1k , j,ik , j,i1k , j,i

z

k ,1 j,ik , j,ik ,1 j,i

y

k , j,1ik , j,ik , j,1i

xk , j,i 

En divisant par le volume de la tranche V = Δx × Δy × Δz, on a une variation de concentration

⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

Δ

+−

ΔΔ+

Δ

+−

ΔΔ+

Δ

+−

ΔΔΔ=Δ +−+−+−

2

1k , j,ik , j,i1k , j,iz

2

k ,1 j,ik , j,ik ,1 j,iy

2

k , j,1ik , j,ik , j,1ixk , j,i

z

CC2C

yx

S

y

CC2C

zx

S

x

CC2C

zy

StK C

 or  yxS;zxS;zyS zyx ΔΔ=ΔΔ=ΔΔ= d’où

⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

Δ

+−+

Δ

+−+

Δ

+−Δ=Δ +−+−+−

2

1k , j,ik , j,i1k , j,i

2

k ,1 j,ik , j,ik ,1 j,i

2

k , j,1ik , j,ik , j,1i

k , j,iz

CC2C

y

CC2C

x

CC2CtK C  

En reprenant les considérations précédentes sur les dérivées, on aboutit alors à :

CK z

C

y

C

x

CK 

t

C2

2

2

2

2

2

Δ=⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂+

∂+

∂=

∂∂

 

La somme des dérivées secondes est appelée Laplacien et notée Δ. On a alors

Page 11: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 11/85

Modélisation

FN CRES

11

CK t

CΔ=

∂∂

 

Il s'agit de l'équation de la diffusion 3D. Ici, il faut bien noter que le coefficient de diffusion

est constant dans toutes les directions de l'espace.Entre le modèle 1D et 3D, on peut bien sûr établir un modèle 2D.

2.4. Conditions limites et condition initialeDans cette démarche, on a travaillé avec une tranche i,j,k "interne", mais il est évident que

  pour la première et la dernière tranche (dans chacune des directions de l'espace), il faudrareprésenter ce qui se passe à chacune des extrémités pour établir un bilan massique du produit

dans la première et dernière tranche.Autrement dit, il faut savoir ce qui se passe sur les frontières du domaine étudié pour pouvoir 

établir un bilan des tranches périphériques. Sans cette connaissance, on ne pourra pas établir la variation de concentration sur ces tranches périphériques, et en cascade, on ne pourra pasétablir de bilan sur une tranche quelconque.

Il s'agit des conditions limites qui gèrent le phénomène sur toutes les frontières spatiales.On distingue généralement 2 grands types de conditions limites :

-  les conditions de type Dirichlet : elles fixent une valeur de la concentration sur la

frontière. Cette condition permet donc toujours de faire un bilan sur la première ou

dernière tranche en établissant la quantité de matière qui franchit la frontière. Dansl'exemple précédent, cela revient à connaitre la concentration sur la frontière du

domaine (par exemple une source de pollution qui maintient constante la

concentration sur la frontière).

-  Les conditions de type Neumann : elles fixent la valeur de la dérivée de la

concentration sur la frontière. Cette condition permet aussi d'établir un bilan sur la

 première ou dernière tranche puisqu'on a directement le terme x/C ∂∂ (ou y/C ∂∂ ou

z/C ∂∂ ). On les appelle aussi condition de flux. Dans l'exemple précédent, cela revient

à avoir du polluant qui entre sans arrêt dans (ou sort de) notre domaine d'étude si ce

flux n'est pas nul).

Insistons sur la différence des 2 conditions par rapport à notre exemple.

Dans le premier cas, du polluant entrera dans (ou sortira de) notre domaine s'il y a unedifférence de concentration entre la frontière et la tranche périphérique (puisque la quantité

qui circule entre 2 tranches est proportionnelle à la différence de concentration). Si lesconcentrations à la frontière et dans la tranche périphérique sont identiques, rien ne rentre

(rien ne sort)

Dans le deuxième cas, quelque soit les concentrations sur la frontière et dans la tranche périphérique, du polluant entre ou sort si le flux est différent de 0, et rien ne se passe si le flux

est nul.

De la même façon, à l'origine temporelle du phénomène, la concentration du produit à

l'intérieur de chaque tranche va commander l'échange de matière entre les tranches, du moins pendant les premiers instants. La connaissance des concentrations à l’instant « zéro » est donc

Page 12: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 12/85

Modélisation

FN CRES

12

indispensable pour appréhender l’évolution des concentrations de chaque tranche dans le

temps. Il s'agit ici de la condition initiale.

Ces conditions sont extrêmement importantes :

1)  elles conditionnent la valeur de la solution, c'est-à-dire qu'elles sont aussi importantesque l'équation elle-même. Des conditions limites différentes, pour une même équation

résolue, donneront des solutions différentes.

2)  elles doivent être compatibles avec l'existence d'une solution, et que celle-ci soitunique (dans le cas contraire, on parle de problème mal posé)

2.5. Régime permanentOn parle d’un phénomène en régime permanent quand son évolution est indépendante du

temps.

Dans l’équation de la diffusion, le termet

C

∂∂ est nul et l'équation devient ΔC = 0 .

Il s’agit de l’équation de Laplace (somme des dérivées secondes en espace nulle).

2.6. Compléments à l’équation de baseOn se remet dans le cas où on regarde ce qui de passe le long de l’axe Ox uniquement (c'est

 plus simple à écrire et l'extension à 2 ou 3 dimensions ne pose pas de problème particulier).

2.6.1. Injection (ou soutirage)

On peut "injecter" ou "soutirer" dans la tranche j une quantité M0 du produit en question.

La quantité totale qui entre ou qui sort de la tranche i s'écrit : 02i1ii MMMM +−= (M0 est

 positif ou négatif selon le cas).

D'oùxS

Mt

x

CC2CK 

xS

M

xS

MM

xS

MC 0

2

1ii1i02i1ii

Δ×+Δ

Δ

+−=

Δ×+

Δ×

−=

Δ×=Δ +−  

SoittxS

M

x

CC2CK 

t

C 0

2

1ii1i

Δ×Δ×+

Δ

+−=

ΔΔ +−  

Ou encore quand Δx et Δt tendent vers zéro : 02

2

q

x

CK 

t

C+

∂=

∂ 

 Remarque : q0 est exprimé en mg/l/s

2.6.2. Cinétique propre du produitIl s'agit d'un terme de croissance (apparition) ou de décroissance (disparition) du produit. On

considère souvent que la variation de concentration due à cette cinétique est proportionnelle àla concentration de produit (cinétique du premier ordre) et se traduit au niveau du bilan de la

tranche i par un terme k'C (en positif ou en négatif).

C'k xS

MMC'k 

xS

MC 2i1ii +

Δ×

−=+

Δ×=Δ d'où

t

C'k 

x

CC2CK 

t

C2

1ii1i

Δ+

Δ+−

=ΔΔ +−  

Page 13: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 13/85

Modélisation

FN CRES

13

Ou encore quand Δx et Δt tendent vers zéro : kCx

CK 

t

C2

2

+∂

∂=

∂∂

 

 Remarque : k est exprimé s-1

 

2.6.3. Convection Nous avons fait l'hypothèse que la matière dans laquelle on observe le produit était immobile.

Supposons maintenant qu'elle est animée d'une vitesse U positive dans le sens des x.

De par ce mouvement, indépendamment du phénomène de diffusion, la quantité de matière

qui entre dans la tranche i en provenance de la tranche i-1 s'écrit :

-  quantité d’eau qui entre par unité de temps : US (on a des m3/s si U en m/s et S en m²)

-  masse de produit par unité de temps : US × Ci-1 (on a des mg/s, si C est en mg/m3)

-  masse de produit : US × Ci-1 × Δt (en mg si Δt en secondes)

et il en sort une quantité USCiΔt par un même raisonnement.

On aura donc le bilan global (accumulation) suivant de matière dans la tranche i :tUSCtUSCMMM i1i2i1ii Δ−Δ+−= −  

et la variation de concentration suivante : tx

CCU

xS

MMC

 j1i2i1i ΔΔ

−+

Δ×

−=Δ −

 

soitx

CCU

x

CC2CK 

t

C i1i

2

1ii1i

Δ−

+−=

ΔΔ −+−  

Le termex

CC i1i

Δ−− est une approximation de la dérivée première en i (au signe près) et tend

donc vers cette dérivée quand Δx → 0. Dans ce cas (Δx→ 0), on a alors l'équation :

x

CU

x

CK 

t

C2

2

∂−

∂=

∂ou bien

2

2

x

CK 

x

CU

t

C

∂=

∂+

∂ 

Il s'agit de l'équation dite de Diffusion-Convection 1D.

En 3D, on considère les vitesses Ux, Uy, Uz, dans chaque direction de l’espace. On établit les

mêmes bilans que précédemment et on obtient :

CK z

CU

y

CU

x

CU

t

Czyx Δ=

∂∂

+∂∂

+∂∂

+∂∂

 

En considérant les vecteurs

-→

U de composantes )U,U,U( zyx  

- Cgrad ⎯→ ⎯ 

de composantes ⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂

∂∂

∂∂

z

C,

y

C,

x

on a l'équation CK Cgrad.Ut

CΔ=+

∂∂ ⎯→ ⎯ →

oùz

CU

y

CU

x

CUCgrad.U zyx ∂

∂+

∂∂

+∂∂

= ⎯→ ⎯ →

 

2.6.4. Cas où K est variable (dans toutes les directions et en tous points)Reprenons le problème sans injection, cinétique ou convection.

Dans un premier temps, on regarde ce qui se passe le long de l'axe Ox.La quantité de matière Mx qui franchit une section S du parallélépipède est :

Page 14: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 14/85

Modélisation

FN CRES

14

tx

CCS)z,y,x(K M 21

x ΔΔ

−=  

tx

CS)z,y,x(K M x Δ

∂∂

−= quand Δx → 0

Introduisons qx le flux, toujours selon l'axe Ox, par unité de surface d'échange :

xxx

x qtSMsoitx

C)z,y,x(K 

tS

Mq Δ=

∂∂

−=Δ

=  

En faisant un bilan d'accumulation dans une tranche i, donc un bilan des flux entre les faces 1

et 2 :

Accumulation = Entrée – Sortie

Accumulation = ( )2x1x qqtS −Δ  

D'où la variation de concentration au sein de la tranche j :

( )2x1x qq

x

t

xS

onAccumulatiC −

Δ

Δ=

Δ×

=Δ  

Où encore

x

qq

x

qq

t

C 1x2x2x1x

Δ

−−=

Δ

−=

ΔΔ

 

En passant en notation différentielle quand Δx → 0 et Δt → 0 :

x

q

t

C x

∂−=

∂∂

avecx

C)z,y,x(K q x ∂

∂−=  

On procède de la même façon sur les axes Oy et Oz.

Le flux sur chacun des axes s'écrit :yC)z,y,x(K q y ∂

∂−= etzC)z,y,x(K q z ∂

∂−=  

En présentation vectorielle, on a

Cgrad).z,y,x(K q −=  

où q est un vecteur de coordonnées )zyx q,q,q

Cgrad ⎯→ ⎯ 

est un vecteur de coordonnées ⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

z

C;

y

C;

x

L'accumulation sur chacun des axes conduit à :

y

q

tC y

∂−=∂∂

etz

qtC z

∂∂−=∂∂ 

La variation totale de concentration sur un volume infinitésimal sera la somme des variations

selon les différents axes :

z

q

y

q

x

q

t

C zyx

∂−

∂−

∂−=

∂∂

 

En notant )q(divr

- lire divergent de q - la somme des dérivées premières :

z

q

y

q

x

q)q(div zyx

∂+

∂+

∂=

Page 15: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 15/85

Modélisation

FN CRES

15

( )Cgrad).z,y,x(K div)q(divt

C−−=−=

∂∂ r

 

C'est l'écriture générale de l'équation de la diffusion (sans convection, ni apport ou soutirage,

ni cinétique)

 Remarque : si K est constant : K(x,y,z) = K, on a alors

( ) ( )CgraddivK CgraddivK )q(divt

C=−−=−=

∂∂ r

 

C.K ²x

²y

²x

C²K 

z

C

zy

C

yx

C

xK 

t

CΔ=⎟⎟

 ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂

+∂∂

+∂∂

=⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ ⎟

 ⎠

 ⎞⎜⎝ 

⎛ ∂∂

∂∂

+⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂

∂∂

+⎟ ⎠

 ⎞⎜⎝ 

⎛ ∂∂

∂∂

=∂∂

 

2.6.5. Autres problématiques Nous avons traité ici le cas d'un produit (dans de l'eau à priori) à travers sa concentration,

mais la démarche peut bien sûr s'appliquer à toute variable : chaleur (c'est d'ailleurs dans cecadre que l'équation de la diffusion fut tout d'abord établie ; on parle alors d'équation de

Fourier) – la variable étudiée est alors la température –, volume d'eau "à surface libre" – la

variable étudiée est alors la hauteur d'eau à travers un bilan des volumes d'eau – …

Dans tous les cas, on arrive à des équations aux dérivées partielles (EDP), généralement

d’ordre 2, c’est-à-dire qu’elles font intervenir les dérivées secondes.

Page 16: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 16/85

Modélisation

FN CRES

16

 

3. CLASSIFICATION DES EQUATIONS AUX DERIVEESPARTIELLES D’ORDRE 2

3.1. Classification par les coniques

La classification des EDP peut être rapprochée de celle des coniques. C’est ce que nous allons

faire.

Soit l’équation d’une conique

ax² + 2bxy + cy² + dx + ey + f = 0 (a ≠ 0 et c ≠ 0)

3.1.1. Classification par discriminant et par valeurs propresLes directions asymptotiques sont obtenues quand on fait tendre les variables x et y vers

l’infini. L’équation d’une conique est alors équivalente aux termes de plus haut degré :

ax² + 2bxy + cy² = 0

En posant y = mx, on obtient

ax² + 2bmx² + cm²x² = 0

cm² + 2bm + a = 0

C’est l’équation des directions asymptotiques. La nature d’une conique dépend du signe du

discriminant Δ = 4b² - 4ac de cette équation :

- si Δ > 0 : 2 directions asymptotiques réelles et distinctes ; la conique est une hyperbole

- si Δ = 0 : 2 directions asymptotiques réelles et confondues ; la conique est une parabole

- si Δ < 0 : pas de directions asymptotiques réelles; la conique est une ellipse

Le signe de Δ est indépendant de tout changement d’axe : c’est une caractéristique intrinsèque

de l’équation.

L’équation des directions asymptotiques ax² + 2bxy + cy² = 0 peut se mettre sous la forme

matricielle suivante :

[ ] 0y

x

c b

 bayx =⎥

⎤⎢⎣

⎡⎥⎦

⎤⎢⎣

⎡soit [ ] [ ] 0yxQyx

T = où ⎥⎦

⎤⎢⎣

⎡=

c b

 baQ

Soit λ1 et λ2 les valeurs propres de Q.

Ce sont les solutions de l’équation (a – λ)(c – λ) – b² = 0

On montre que la classification des coniques par discriminant est équivalente à :

- si λ1 et λ2 sont non nulles et de signes différents ; la conique est une hyperbole

- si λ1 ou λ2 sont nulles ; la conique est une parabole

- si λ1 et λ2 sont non nulles et de même signe ; la conique est une ellipse

3.1.2. Application aux EDP d’ordre 2Les propos précédents peuvent s’étendre directement aux équations aux dérivées partielles

linéaires d'ordre 2.

Page 17: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 17/85

Modélisation

FN CRES

17

Soit l’EDP linéaire générale à 2 variables :

0)y

u,

x

u,u,y,x(F

y

uc

yx

u b

x

ua

2

22

2

2

=∂∂

∂∂

+∂∂

+∂∂

∂+

∂∂

où a, b et c sont des fonctions de x et y

On considère la matrice ⎥⎦⎤⎢

⎣⎡= c b

 baQ et ses valeurs propres.

L’EDP est hyperbolique si les valeurs propres sont non nulles et de signes différents.

Par exemple, l’équation des ondes2

2

2

2

2

2

2

22

t

u

z

u

y

u

x

uk 

∂∂

=⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂

+∂∂

+∂∂

 

L’EDP est parabolique si au moins une valeur propre est nulle.

Par exemple, l’équation de la diffusiont

u

z

u

y

u

x

uk 

2

2

2

2

2

2

∂∂=⎟⎟

 ⎠

 ⎞⎜⎜⎝ 

⎛ ∂∂+

∂∂+

∂∂  

L’EDP est elliptique si les valeurs propres sont non nulles et de même signe.

Par exemple, l’équation de Laplace 0z

u

y

u

x

u2

2

2

2

2

2

=∂∂

+∂∂

+∂∂

 

 Remarque : dans la mesure où a, b et c dépendent de x et y, l’EDP peut être de type mixte.

3.2. Equation stationnaire et équation d’évolutionOn a vu que les équations hyperboliques et paraboliques avaient des directions asymptotiques

réelles. Cela veut dire que l’une des variables peut tendre vers l’infini. Aussi, toute EDP

  physique hyperbolique ou parabolique fait intervenir le temps. On parle alors d’équations

évolutives, car la solution de ces équations évolue en fonction du temps.

A contrario, la solution d’une équation elliptique n’ayant pas de directions asymptotiques

réelles, le temps ne peut pas y intervenir. On parle d’équations stationnaires.

Page 18: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 18/85

Modélisation

FN CRES

18

 

4. RESOLUTION DES EDP - PRINCIPE DES METHODESNUMERIQUES

On aborde à partir de ce chapitre les méthodes de résolution des EDP. Celles-ci n’ont

généralement pas de solutions exactes et on utilise donc des méthodes numériques.

Les méthodes numériques s’intéressent à la recherche de valeurs de la fonction en des

endroits particuliers. Autrement dit, on ne cherche pas l’écriture d’une fonction qui vérifie

l’équation, mais par quelles valeurs passe la fonction en des abscisses particulières (c’est la

méthode des différences finies), ou bien on recherche sur des éléments du domaine étudié

l’écriture d’une fonction simple qui approxime au mieux la solution recherchée (c’est la

méthode des éléments finis).

Une fois ce travail fait, on a donc une « image graphique » de la solution.

 Nous allons étudier successivement ces 2 méthodes.

Page 19: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 19/85

Modélisation

FN CRES

19

 

5. LA METHODE DES DIFFERENCES FINIES

5.1. Principe de la méthode

Ce principe se décline en plusieurs étapes.

a - Le domaine étudié est maillé. On parle de discrétisation du domaine.

Exemple de maillage d’un domaine

1D. Ce maillage peut être régulier 

ou non, c’est-à-dire que le pas du

maillage peut être constant ou non

Frontière

 NoeudsPas du maillage

Exemple de maillage d’un

domaine 2D (on n’en a représenté

qu’une partie)

Frontière

Maillage

 Noeuds

La méthode des différences finies recherche une solution aux nœuds du maillage.

  b – On discrétise également l’EDP, c’est-à-dire qu’on va écrire en chaque nœud une

approximation algébrique de l’équation d’origine.

c – On écrit autant d’équations algébriques qu’il y a de nœuds où on cherche une solution, ce

qui conduit à écrire un système d’équations 

d – on résout ce système d’équations 

 Nous traiterons différents exemples simples illustratifs de cette méthode.

Page 20: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 20/85

Modélisation

FN CRES

20

5.2. Cas d’une EDP elliptique (stationnaire)

5.2.1. Cas de conditions de Dirichlet

On considère un domaine carré, donc 2

dimensions d’espace (problème 2D), oùon veut résoudre le problème suivant :

0y

u

x

u2

2

2

2

=∂

∂+

∂ 

(équation de Laplace)

Les conditions limites sont de type

Dirichlet, telles que décrites sur le schéma

ci-contre, g1 et g2 étant 2 constantes

quelconques.

0y

u

x

u2

2

2

2

=∂

∂+

∂ 

u = g2 

u = g1 u = 0

u = 0 x

y

 

a – Discrétisation du domaine d’étude

Le domaine d’étude comprend le domaine interne, noté Ω, et la frontière, notée Γ.

On va créer un maillage relativement lâche (afin de limiter ensuite les écritures ; mais

l’extension à plus de nœuds ne pose aucun problème) et régulier.

Ce maillage introduit en tout 25 nœuds,dont 9 nœuds internes et 16 nœuds

frontières.

Cependant, les nœuds frontières ne

constituent pas des valeurs à rechercher,

  puisque la valeur de la fonction y est

connue (conditions de Dirichlet).

Seuls les nœuds ici numérotés de 1 à 9

constituent donc les inconnues de notre

  problème, sachant que bien sûr, les

conditions limites doivent intervenir dansla solution.

u = g2 

u = g1 u = 0

u = 0 x

y

1 2 3

4 5 6

7 8 9

Ω

Γ 

Le maillage étant régulier, on pose Δx = Δy = h

Rappelons que la méthode des différences finies permet de trouver une valeur approchée de la

solution aux 9 nœuds, solution qui vérifie l’équation 0y

u

x

u2

2

2

2

=∂

∂+

∂ 

b – Discrétisation de l’équationCette étape consiste à remplacer l’EDP par une équation algébrique approchée. Plusieurs

approches sont possibles.

Page 21: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 21/85

Modélisation

FN CRES

21

Approche graphique (plus pédagogique)

On a déjà vu qu’on pouvait remplacer une dérivée par une équation. Considérons la courbe ci-

dessous et sa dérivéex

u

∂∂

au point d'abscisse x. Cette dérivée correspond à la tangente à la

courbe au point d'abscisse x.

On note x = xi , x + Δx = xi + 1 et x - Δx = xi - 1 u(x) = ui , u(x + Δx) = ui + 1 et u(x - Δx) = ui - 1 

On peut approximer cette tangente

  par la pente de la corde entre les

 points (xi ,ui) et (xi+1 ,ui+1).

Cette pente est égale àx

uu i1i

Δ

−+ ;

 plus Δx est petit, et plusl'approximation est valide.

Autrement dit, le termex

uu i1i

Δ

−+ est

une approximation de la dérivée

 premièreix

u⎟

 ⎠

 ⎞⎜⎝ 

⎛ ∂∂

à l'abscisse xi quand Δx→ 0.

Cette discrétisation de la dérivée premièrex

uu

x

u i1i

Δ

−≈

∂∂ + est dite avant, ou à droite, ou

 progressive, ou aval, car elle fait intervenir ui+1

à l’abscisse xi+1

.

De la même façon, on peut approximer la dérivée premièreix

u⎟

 ⎠

 ⎞⎜⎝ 

⎛ ∂∂

à l'abscisse xi par la pente

de la corde entre les points (xi ,ui) et (xi-1 ,ui-1) :x

uu

x

u 1ii

i Δ

−≈⎟

 ⎠

 ⎞⎜⎝ 

⎛ ∂∂ − quand Δx→ 0.

Cette discrétisation de la dérivée premièrex

uu

x

u 1ii

Δ

−≈

∂∂ − est dite arrière, ou à gauche, ou

régressive, ou amont, car elle fait intervenir ui-1 à l’abscisse xi-1.

On peut également approximer la dérivée première par la pente de la corde entre les abscisses

xi-1 et xi+1. On a alorsx2

uu

x

u 1i1i

Δ

−≈

∂∂ −+  

Cette discrétisation de la dérivée premièrex2

uu

x

u 1i1i

Δ

−≈

∂∂ −+ est dite centrée.

On peut maintenant avoir une approximation de la dérivée seconde en x i , en dérivant 2 fois

  par rapport à x, et en appliquant une fois une discrétisation avant, puis une discrétisation

arrière :

ui

x

xi xi +1xi -1

Δx Δx

ui +1

u

tangente

Approximationde la tangente

ui -1

Page 22: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 22/85

Modélisation

FN CRES

22

2

1ii1i

2

1ii

2

i1i

i1i

i1i

ii

2

2

x

uu2u

x

uu

x

uu

x

u

xx

u

xx

uu

xx

u

xx

u

Δ

+−=

Δ

−−

Δ

−=⎟

 ⎠

 ⎞⎜⎝ 

⎛ Δ∂

∂−⎟

 ⎠

 ⎞⎜⎝ 

⎛ Δ∂

∂=⎟

 ⎠

 ⎞⎜⎝ 

⎛ Δ

∂∂

≈⎟ ⎠

 ⎞⎜⎝ 

⎛ ∂∂

∂∂

=⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂ +−−+

+

+

 

Cette discrétisation de la dérivée seconde2

1ii1i

2

2

x

uu2u

x

u

Δ

+−≈

∂ +− est dite centrée.

Approche analytique (plus rigoureuse)

On peut arriver au même résultat en utilisant le développement de Taylor d’une fonction au

voisinage d’un point.

...x

)x(u

!4

x

x

)x(u

!3

x

x

)x(u

!2

x

x

)x(ux)x(u)xx(u

4

44

3

33

2

22

+∂

∂Δ+

∂∂Δ

+∂

∂Δ+

∂∂

Δ+=Δ+  

ou encore ...x

)x(u

!4

x

x

)x(u

!3

x

x

)x(u

!2

x

x

)x(ux)x(u)xx(u

4

44

3

33

2

22

+∂

∂Δ+

∂∂Δ

−∂

∂Δ+

∂∂

Δ−=Δ−  

 Remarque : les 3 points de suspension à la fin des expressions précédentes représentent une

quantité que l'on néglige dans l'égalité quand Δ x tend vers zéro. C'est-à-dire que ce terme

« tend plus rapidement » vers zéro que le dernier terme de la série écrit.

Cette quantité négligée est écrite sous forme d'une fonction o(Δ xn) où n est alors appelé 

"ordre de développement". o(Δ xn) se lit alors comme un "terme de l'ordre de Δ x

n".

 Ex : )x(ox

)x(u

!4

x

x

)x(u

!3

x

x

)x(u

!2

x

x

)x(ux)x(u)xx(u 5

4

44

3

33

2

22

Δ+∂

∂Δ+

∂Δ+

∂Δ+

∂∂

Δ+=Δ+  

Plus l'ordre du terme négligé est élevé, et plus l'égalité écrite est "vrai", c'est-à-dire que plus

le terme négligé est petit; donc mieux c'est !

En écrivant ces développements à l’ordre 1 :

)x(ox

)x(ux)x(u)xx(u 2Δ+

∂∂

Δ+=Δ+  

soit )x(ox

)x(u)xx(u

x

)x(uΔ+

Δ−Δ+

=∂

∂; c’est le schéma de discrétisation avant.

)x(ox

)x(ux)x(u)xx(u 2Δ+∂

∂Δ−=Δ−  

soit )x(ox

)xx(u)x(u

x

)x(uΔ+

ΔΔ−−

=∂

∂; c’est le schéma de discrétisation arrière.

Pour le schéma centré :

)x(ox

)x(u

!2

x

x

)x(ux)x(u)xx(u 3

2

22

Δ+∂

∂Δ+

∂∂

Δ+=Δ+  

)x(ox

)x(u

!2

x

x

)x(ux)x(u)xx(u

3

2

22

Δ+∂

∂Δ+∂

∂Δ−=Δ−  

Page 23: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 23/85

Modélisation

FN CRES

23

en faisant la différence de ces 2 expressions :

)x(ox

)x(ux2)xx(u)xx(u 3Δ+

∂∂

Δ=Δ−−Δ+  

soit )x(ox2

)xx(u)xx(u

x

)x(u 2Δ+Δ

Δ−−Δ+=

∂∂

; schéma de discrétisation centré

 Remarque : le schéma centré est meilleur que les 2 autres, puisque avec une erreur en o(Δ x²)

alors que les schémas avant et arrière ont une erreur en o(Δ x).Ceci se visualise aisément par 

la méthode graphique reprise ci-dessous.

On "voit" que quand Δ x tend vers zéro, l'approximation de la tangente en schéma centré tend 

 plus rapidement vers la bonne tangente qu'en schéma décentré.

ui

x

xi xi +1

Δx

ui +1

u

tangente

Approximationde la tangente

ui

x

xi xi +1xi -1

Δx Δx

ui +1

u

tangente

Approximationde la tangente

ui -1

Schéma de discrétisation avant :

erreur en o(Δ x)

Schéma de discrétisation centré :

erreur en o(Δ x²)

Pour la dérivée seconde :

)x(ox

)x(u

!3

x

x

)x(u

!2

x

x

)x(ux)x(u)xx(u 4

3

33

2

22

Δ+∂

∂Δ+

∂Δ+

∂∂

Δ+=Δ+  

)x(ox

)x(u

!3

x

x

)x(u

!2

x

x

)x(ux)x(u)xx(u 4

3

33

2

22

Δ+∂

∂Δ−

∂Δ+

∂∂

Δ−=Δ−  

en faisant la somme de ces 2 expressions :

)x(ox

)x(ux)x(u2)xx(u)xx(u 4

2

22 Δ+

∂Δ+=Δ−+Δ+  

soit )x(ox

)xx(u)x(u2)xx(ux

)x(u 222

2

Δ+Δ Δ−+−Δ+=∂∂ ; schéma centré

En revenant à notre exemple et en notant

u(x,y) = ui , j , u(x ± Δx , y) = ui ± 1 , j , u(x, y ± Δy) = ui , j ± 1 , u(x ± Δx , y ± Δy) = ui ± 1 , j ± 1 

Par ailleurs, nous rappelons que le maillage est régulier : Δx = Δy = h

)h(o

h

uu2u

x

u 2

2

 j,1i j,i j,1i

2

2

++−

=

∂ +− 

Page 24: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 24/85

Modélisation

FN CRES

24

)h(oh

uu2u

y

u 2

2

1 j,i j,i1 j,i

2

2

++−

=∂

∂ +− 

En sommant ces 2 expressions

)h(o

h

uu2u

h

uu2u

y

u

x

u 2

2

1 j,i j,i1 j,i

2

 j,1i j,i j,1i

2

2

2

2

++−

++−

=

∂+

∂ +−+− 

)h(oh

uuu4uu

y

u

x

u 2

2

1 j,i j,1i j,i1 j,i j,1i

2

2

2

2

+++−+

=∂

∂+

∂ ++−− 

Comme l'équation à résoudre est 0y

u

x

u2

2

2

2

=∂

∂+

∂ 

0h

uuu4uu2

1 j,i j,1i j,i1 j,i j,1i =++−+ ++−−

en négligeant l’écriture de o(h²)

soit 0uuu4uu 1 j,i j,1i j,i1 j,i j,1i =++−+ ++−−

 Cette expression constitue donc une approximation de l’EDP écrite au nœud ui,j . C’est une

relation entre 5 nœuds qui approxime l’EDP écrite au nœud u i,j. Calquée sur le maillage, on

l’écrit sous forme graphique comme ci-dessous à droite où on représente les coefficients de

chacun des nœuds :

Ui,j

Ui,j-1

Ui,j+1

Ui+1,jUi-1,j-4 11

1

1

 Discrétisation de l'équation de Laplace

c – Ecriture du système d’équationsOn va maintenant écrire l’équation discrétisée en chacun des nœuds inconnus, ce qui va

conduire à un système de 9 équations à 9 inconnues. Pour cela, on reprend la numérotation

des 9 nœuds de 1 à 9 ; en effet, cela est moins laborieux que de travailler en double

coordonnées (i,j).

Page 25: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 25/85

Modélisation

FN CRES

25

  u = g2 

u = g1u = 0

u = 0 x

y

1 2 3

4 5 6

7 8 9

nœud 1 : 0uuu4 421 =++−  

les nœuds à gauche et en bas correspondent

aux conditions de Dirichlet nulles

nœud 2 : 0uuu4u 5321 =++−  

nœud 3 : 0ugu4u 6132 =++−  le nœud à droite correspond à la condition de

Dirichlet égale à g1

nœud 4 : 0uuu4u 7541 =++−  

nœud 5 : 0uuu4uu 86524 =++−+  

nœud 6 : 0ugu4uu 91635 =++−+  

nœud 7 : 0guu4u 2874 =++−  

nœud 8 : 0guu4uu 29857 =++−+  

nœud 9 : 0ggu4uu 21968 =++−+  

On représente généralement ce système d’équations sous forme matricielle, ce qui permet de

 bien choisir ensuite la méthode de résolution qui sera adoptée. Ici, la première ligne sert à

  bien repérer les colonnes, sachant que l’absence d’une valeur dans la matrice carrée

correspond à zéro.

⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎢⎢⎢

−−

=

⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎢⎢⎢

⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎢⎢⎢

21

2

2

1

1

9

8

7

6

5

4

3

2

1

987654321

gg

g

g

g

0

0

g

0

0

u

u

u

u

u

u

u

u

u

411

1411

141

1411

11411

1141

141

1141

114

uuuuuuuuu

 

d – Résolution du système d’équations Nous renvoyons le lecteur aux méthodes numériques de résolution d’un tel système (dernier 

chapitre). Nous l’avons résolu avec Excel dans l’application numérique suivante : le carré fait

1 mètre de coté (donc h=0.25m) et avec g1 = 1 , g2 = 2.

 Nous obtenons l’image suivante de la solution :

Page 26: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 26/85

Modélisation

FN CRES

26

  u = 2

u = 1u = 0

u = 0 x

y

0.15 0.34 0.55

0.25 0.67 0.87

0.86 1.20 1.27

Un maillage plus fin permettrait d’obtenir une meilleure précision (mais avec plus d’équations

à résoudre et au détriment de l’approche pédagogique).

5.2.2. Cas de conditions de Neumann homogènes

On considère le même problème que

 précédemment : un domaine carré, donc 2

dimensions d’espace, où on veut résoudre

le problème suivant : 0y

u

x

u2

2

2

2

=∂

∂+

∂ 

Les conditions limites sont de type

Dirichlet à droite et à gauche, telles que

décrites sur le schéma ci-contre, et de type

 Neumann homogène (c'est-à-dire nulle) en

haut et en bas.

0y

u

x

u2

2

2

2

=∂

∂+

∂ 

0y

u=

∂∂

u = g2 u = g1

0y

u=

∂∂ x

y

 

a – Discrétisation du domaine d’étude

Page 27: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 27/85

Modélisation

FN CRES

27

 

Afin de limiter le nombre d’équations à

écrire, on adopte le maillage suivant du

dessin de droite.

Les nœuds inconnus y sont numérotés de 1à 8. S’il n’y a pas de problème pour les

nœuds internes 3 à 6, on peut

légitimement s’interroger sur les nœuds

frontières 1, 2 et 7, 8. En effet, la valeur de

la dérivée est connue en ces nœuds, mais

en aucun cas la valeur de la fonction. Ce

sont donc bien des nœuds inconnus.

0y

u=

∂∂

u = g2 u = g1

0y

u=

∂∂ x

y

1 2

3 4

5 6

7 8

b – Discrétisation de l’équationOn a exactement la même équation que précédemment ; elle peut être approximée par 

l’expression 0uuu4uu 1 j,i j,1i j,i1 j,i j,1i =++−+ ++−− écrite au nœud inconnu ui,j

Ceci introduit une difficulté aux nœuds frontières puisque cette discrétisation de l’équation

introduit à chaque fois un nœud en dehors du domaine d’étude, qu’on appelle « nœud fictif ».

Ui,j

Ui,j-1

Ui,j+1

Ui+1,jUi-1,j

Frontière inférieure

 Nœud fictif 

Ω  

Ui,j

Ui,j-1

Ui,j+1

Ui+1,jUi-1,j

Frontière supérieure

 Nœud fictif 

Ω  

Cas de la frontière inférieure Cas de la frontière supérieure

Pour éliminer les nœuds fictifs de la discrétisation, il suffit de prendre en compte la condition

de flux sur la frontière en question.

Prenons l'exemple du nœud n°1 et nommons le nœud fictif 3'.•  La discrétisation de l'équation au nœud 1 donne :

0uuu4ug 321'31 =++−+  

•  La discrétisation centrée de la condition de flux au nœud 1 s'écrit :

0h2

uu

y

u '33 =−

≈∂∂

ce qui implique que u3' = u3 

On remplace u3' dans la discrétisation de l'équation

0uuu4ug 32131 =++−+ soit 0u2uu4g 3211 =++−  

Ceci qui permet d'éliminer le nœud fictif de l'équation écrite au nœud 1.

Page 28: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 28/85

Modélisation

FN CRES

28

Remarquons que cette discrétisation vérifie simultanément l'équation et la condition limite.

On obtient donc 3 formes différentes de discrétisation selon la position du nœud :

 Nœuds de la frontière

supérieure

-4 11

2

 Nœuds internes

-4 11

1

 Nœuds de la frontière

inférieure

-4 11

2

c – Ecriture du système d’équations

0y

u=

∂∂

u = g2 u = g1

0y

u=

∂∂ x

y

1 2

3 4

5 6

7 8  Nœud 1 : 0u2uu4g 3211 =++−  

 Nœud 2 : 0u2gu4u 4221 =++−  

  Nœud 3 : 0uuu4ug 54311 =++−+  

  Nœud 4 : 0ugu4uu 62423 =++−+  

  Nœud 5 : 0uuu4ug 76531 =++−+  

  Nœud 6 : 0ugu4uu 82645 =++−+  

  Nœud 7 : 0uu4u2g8751

=+−+  

  Nœud 8 : 0gu4u2u 2867 =+−+  

La mise sous forme matricielle donne :

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

−−

=

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

−−

2

1

2

1

2

1

2

1

8

7

6

5

4

3

2

1

87654321

g

g

g

g

g

gg

g

u

u

u

u

u

uu

u

412

142

1411

1141

1411

1141241

214

uuuuuuuu

 

d – Résolution du système d’équations

Page 29: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 29/85

Modélisation

FN CRES

29

En prenant g1 = 1 et g2 = 2 dans un domaine carré de 1 mètre de côté, nous obtenons l'image

suivante de la solution :

u = 2u = 1 

x

y

1.33

1.33

1.33

1.33

1.67

1.67

1.67

1.67

La solution dans cet exemple est un plan entre les valeurs 1 et 2 selon Ox. Ce qui vérifie parfaitement l'équation à résoudre (toutes les dérivées secondes sont nulles, donc leur somme

est nulle) et les dérivées premières en y sont nulles.

5.2.3. Cas de conditions de Neumann non homogènesIl n'y a pas de difficulté particulière, il suffit de suivre la méthode en intégrant la condition de

flux non homogène au niveau de la discrétisation.

Prenons comme exemple des conditions de flux égales à f 1 à la frontière supérieure et f 2 à la

frontière inférieure.

1f yu =

∂∂

u = g2 u = g1

2f y

u=

∂∂

 x

y

1 2

3 4

5 6

7 8

a – Discrétisation du domaine d’étudeCette étape est strictement identique à la précédente

b – Discrétisation de l’équationOn a exactement la même équation que précédemment ; elle peut être approximée par 

l’expression 0uuu4uu 1 j,i j,1i j,i1 j,i j,1i =++−+ ++−− écrite au nœud inconnu ui,j

Page 30: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 30/85

Modélisation

FN CRES

30

Comme auparavant, des nœuds fictifs apparaissent; on discrétise donc la condition de flux.

Prenons l'exemple du nœud n°1 et nommons le nœud fictif 3'.

•  La discrétisation de l'équation au nœud 1 donne :

0uuu4ug 321'31 =++−+  

•  La discrétisation centrée de la condition de flux au nœud 1 s'écrit :

2'33 f 

h2

uu

y

u=

−≈

∂∂

ce qui implique que 23'3 hf 2uu −=  

On remplace u3' dans la discrétisation de l'équation

0uuu4hf 2ug 321231 =++−−+ soit 0u2uu4hf 2g 32121 =++−−  

c – Ecriture du système d’équations

1f y

u=

∂∂

u = g2 u = g1

2f y

u=

∂∂

 x

y

1 2

3 4

5 6

7 8 Nœud 1 : 23211 hf 2u2uu4g =++−  

 Nœud 2 : 24221 hf 2u2gu4u =++−  

  Nœud 3 : 0uuu4ug 54311 =++−+  

  Nœud 4 : 0ugu4uu 62423 =++−+  

  Nœud 5 : 0uuu4ug 76531 =++−+  

  Nœud 6 : 0ugu4uu 82645 =++−+  

 Nœud 7 : 18751 hf 2uu4u2g −=+−+   Nœud 8 : 12867 hf 2gu4u2u −=+−+  

La mise sous forme matricielle donne :

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

−−

−−

−−

=

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

−−

21

11

2

1

2

1

22

12

8

7

6

5

4

3

2

1

87654321

ghf 2

ghf 2

g

g

g

g

ghf 2

ghf 2

u

u

u

u

u

u

u

u

412

142

1411

1141

1411

1141

241

214

uuuuuuuu

 

Page 31: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 31/85

Modélisation

FN CRES

31

 

5.3. Notions sur les erreurs de la méthode

5.3.1. Présentation

Soit une équation aux dérivées partielles en u notée L(u) = g où L est un opérateur 

différentiel. Par exemple, dans l’équation de Laplace, L est l’opérateur 2

2

2

2

yx ∂∂

+∂∂

, et L(u) = 0

est alors l’équation 0y

u

x

u2

2

2

2

=∂∂

+∂∂

 

On peut alors formaliser l’équation discrétisée sous la forme LΔ(u) = g où Δ est

symboliquement le pas de discrétisation et LΔ un opérateur algébrique.

La solution obtenue par la méthode des différences finies peut alors se noter uΔ , c’est-à-dire

l’approximation de u avec le pas de discrétisation Δ.

LΔ(u)L(u)

uΔ u

Discrétisation

RésolutionCe que l’on cherche

Représentativité

 

On peut étudier les erreurs au niveau des différentes étapes de ce schéma (de A. Le Pouriet)

a)  La discrétisation est liée à la notion de consistance (ou de cohérence). Un schéma de

discrétisation consistant est un schéma qui représente bien l’équation d’origine.

 b)  La résolution des équations obtenues est liée à la notion de stabilité de la solution.

Une solution uΔ stable doit vérifier le problème LΔ(u) = g

c)  La représentativité de uΔ est liée à la notion de convergence. uΔ doit être une

approximation convenable de u.

d) 

LΔ(u)L(u)

uΔ u

 Discrétisation

 Résolution

 Représentativité 

Consistance

Stabilité

Convergence

Page 32: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 32/85

Modélisation

FN CRES

32

 

5.3.2. La consistanceOn appelle « erreur de troncature » la quantité

R Δ(u) = LΔ(u) – L(u)

c’est-à-dire concrètement la différence entre le schéma discrétisé et l’équation d’origine.

LΔ(u) est dit consistant si R Δ(u) 0 quand Δ 0

Reprenons l’exemple précédent en posant L(u) =2

2

2

2

y

u

x

u

∂∂

+∂∂

et étudions la consistance.

L’opérateur différentiel "dérivée seconde par rapport à une variable z" peut se symboliser par 

2

2

Zz

L∂∂

= et l’équation peut alors se représenter par :

2

2

2

2

YXyu

xu)u(L)u(L)u(L

∂∂+∂∂=+=  

Le développement de Taylor au voisinage de x s’écrit :

...x

u

!4

x

x

u

!3

x

x

u

!2

x

x

ux)y,x(u)y,xx(u

4

44

3

33

2

22

+∂∂Δ

+∂∂Δ

±∂∂Δ

+∂∂

Δ±=Δ±  

On en déduit :

...x

u

12

x

x

ux)y,x(u2)y,xx(u)y,xx(u

4

44

2

22 +

∂∂Δ

+∂∂

Δ+=Δ−+Δ+  

...x

u

12

x

x

u

x

)y,xx(u)y,x(u2)y,xx(u4

42

2

2

2 +∂∂Δ

=∂∂

−ΔΔ−+−Δ+

 

4 4 34 4 213214 4 4 4 4 4 4 34 4 4 4 4 4 4 21)u(R )u(L)u(L xxx ΔΔ =−  

De la même façon, le développement de Taylor au voisinage de y s’écrit :

...y

u

!4

y

y

u

!3

y

y

u

!2

y

y

uy)y,x(u)yy,x(u

4

44

3

33

2

22

+∂∂Δ

+∂∂Δ

±∂∂Δ

+∂∂

Δ±=Δ±  

On aboutit à :

...y

u

12

y

y

u

y

)yy,x(u)y,x(u2)yy,x(u

4

42

2

2

2

+∂

∂Δ=

∂−

Δ

Δ−+−Δ+ 

4 4 34 4 213214 4 4 4 4 4 4 34 4 4 4 4 4 4 21)u(R )u(L)u(L yyy ΔΔ =−  

On peut sommer les 2 expressions du dessus :

)u(R )u(R )u(R )u(L)u(L)u(L)u(L)u(L)u(L yxyxyx ΔΔΔΔΔΔ =+=−=−−+  

L’erreur de troncature totale est donc ...y

u

12

y

x

u

12

x)u(R 

4

42

4

42

+∂∂Δ

+∂∂Δ

=Δ  

Cette quantité tend vers 0 quand Δx et Δy tendent vers 0. Le schéma de discrétisation est donc

consistant.

Page 33: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 33/85

Modélisation

FN CRES

33

 

5.3.3. La stabilitéCette notion concerne les résolutions des systèmes d’équations pour lesquels on est souvent

amené à mettre en œuvre des algorithmes itératifs. La stabilité peut se résumer à dire qu’une

erreur d’arrondi (ou une perturbation numérique) ne doit pas s’amplifier au cours des calculs.

On peut imager cette notion de stabilité avec une bille et un bol.On retourne le bol et on pose la bille en équilibre. La moindre

 perturbation de la bille va entraîner sa chute. C’est un système instable

(la perturbation s’amplifie)

A contrario, on met le bol à l’endroit et la bille au fond. Une

 perturbation de la bille va entraîner une oscillation de celle-ci au fond

du bol, mais elle va revenir à un état d’équilibre. C’est un système

stable.

L’étude mathématique de la stabilité est assez complexe et elle n'est pas développée ici. Un

exemple de résolution est proposé dans les exercices illustrant l'instabilité d'un schéma

numérique.

5.3.4. La convergenceLe théorème de Lax dit que s'il y a consistance et stabilité, alors il y a convergence.

Autrement dit, un schéma de discrétisation consistant et une méthode de résolution stable

conduisent à une solution qui est une bonne image de la solution recherchée.

Dans la méthode des différences finies, les discrétisations vues jusqu'à maintenant (les plus

classiques) sont consistantes. Aussi, on s'attachera à préciser uniquement les conditions de

stabilité pour vérifier la convergence.

5.4. Cas d’une EDP parabolique (évolutive)

5.4.1. Position du problème et méthode

 Nous allons étudier l’exemple de l’équation de la diffusion 1D :2

2

x

ua

t

u

∂∂

=∂∂

 

La différence fondamentale par rapport au problème précédent est l’introduction du temps

dans l’équation.

Supposons qu’on travaille entre 0 et L en Ox,

on obtient un domaine non borné en temps si

on cherche à dessiner le domaine où on doit

résoudre l’équation.

0 L

xt0 

t

Ω 

En effet, on commence à étudier le phénomène à t = t0 , mais il n’y a pas de frontière pour t > 0 (en tous cas, on ne peut pas écrire de condition limite à un temps t 1 > t0 car cela n’a pas

Page 34: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 34/85

Modélisation

FN CRES

34

de sens physique : cela voudrait dire que quelque chose qui se passe dans l’avenir influence

quelque chose qui se situe dans le passé – c’est de la science fiction ! – )

Dans un tel cas, le principe de la méthode

consiste, connaissant la solution à un pas de

temps ti à calculer la solution au pas de tempsti+1 , puis connaissant cette solution à ti+1, on

calcule celle à ti+2 et ainsi de suite jusqu’au

moment où on décide d’arrêter les calculs.

On procède donc pas de temps par pas de

temps, la première itération se faisant à partir 

de la condition initiale (à t0). 0 L

x

t0

t

Ω

 

ti 

ti+1 Calcul de la solution à ti+1

ti+2 

Calcul de la solution à ti+2

5.4.2. Les schémas typesOn va considérer des conditions limites de type Dirichlet (des conditions de flux ne posant

 pas de problème particulier en introduisant des nœuds fictifs).

On pose la condition limite à gauche = g1 et celle de droite = g2.

On  rappelle qu’on a aussi une condition initiale à t = t0 et sans ces conditions limites et

initiale, le problème physique n’a pas de solution.

5.4.2.1   Discrétisation du domaine d’étude

Le domaine est discrétisé de façon

classique en espace et en temps. Nousavons représenté ci-contre une partie du

maillage qui peut être régulier ou non en

temps et en espace ; pour simplifier les

écritures, nous le prendrons régulier en

temps et régulier en espace.

0 L

xt0

t

x x+Δxx-Δx

t

t+Δt

u = g1 u = g2

5.4.2.2   Discrétisation de l’équation Nous allons considérer 3 cas pour discrétiser l’équation.

Cas 1

22

2

x

)t,xx(u)t,x(u2)t,xx(u

x

u

ΔΔ++−Δ−

=∂∂

 

t

)t(u)tt,x(u

t

u

Δ

−Δ+=

∂ 

La somme de ces 2 termes de l’EDP à partir de

ces 2 discrétisations introduit une équation qui

relie 4 nœuds du maillage. Cette équation

 permet de calculer ce qui se passe à l’instant t +

Δt en fonction de l’instant t de façon explicite. 0 L

xt0

t

xx+Δ

xx-Δ

x

t

t+Δt

Schéma explicite

Page 35: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 35/85

Modélisation

FN CRES

35

 

En écrivant cette équation aux n nœuds inconnus du temps t + Δt, on obtient ce qui se passe à

cet instant de façon simple (n équations, chacune à 1 inconnue). On parle d’un schéma

explicite car on peut connaître la solution en 1 nœud à t + Δt indépendamment des autres

nœuds à t + Δt.

La condition initiale fournira la première ligne (à t = t 0) pour pouvoir démarrer les calculs.

Cas 2

22

2

x

)tt,xx(u)tt,x(u2)tt,xx(u

x

u

ΔΔ+Δ++Δ+−Δ+Δ−

=∂∂

(il n’y a aucun problème à écrire la

discrétisation de la dérivée en espace à l’instant t + Δt)

t

)t(u)tt,x(u

t

u

Δ

−Δ+=

∂ 

La somme de ces 2 termes de l’EDP à partir deces 2 discrétisations introduit une équation qui

relie toujours 4 nœuds du maillage. Cette

équation ne permet pas de calculer directement

comme précédemment ce qui se passe à l’instant

t + Δt en fonction de l’instant t puisqu’elle relie

3 nœuds inconnues.

0 L

xt0

t

x x+Δxx-Δx

t

t+Δt

Schéma implicite

Il faut donc écrire le système de n équations, chacune à 3 inconnues, et le résoudre, ce qui est

 plus complexe que le schéma explicite. On parle ici d’un schéma implicite dans la mesure où

la solution n’est connue qu’à travers ce système d’équations (on ne peut pas calculer lasolution en 1 nœud sans connaître toute la solution au temps t + Δt).

La condition initiale fournira la première ligne (à t = t 0) pour pouvoir démarrer les calculs.

Cas 3

22

2

x

)t,xx(u)t,x(u2)t,xx(u

x

u

ΔΔ++−Δ−

=∂∂

 

t2

)tt,x(u)tt,x(u

t

u

Δ

Δ−−Δ+=

∂ 

La somme de ces 2 termes de l’EDP à partir de

ces 2 discrétisations introduit une équation qui

relie 5 nœuds du maillage. Cette équation

  permet de calculer ce qui se passe à l’instant

t + Δt en fonction de l’instant t et t - Δt de façon

explicite.

0 L

xt0

t

x x+Δxx-Δx

t

t+Δt

Schéma explicite à 2 pas

t-Δt

On parle d’un schéma explicite à 2 pas. Cependant, elle nécessite de connaître 2 pas de temps

successifs pour démarrer les calculs, or seul un pas de temps est connu avec la condition

initiale. Il faut donc faire une hypothèse sur le premier pas de temps t1, ou bien démarrer les

calculs par un des deux schémas précédents pour connaître les instants t0 et t1.

Page 36: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 36/85

Modélisation

FN CRES

36

 Nous allons maintenant détailler ces 3 schémas.

5.4.3. Le schéma expliciteEcrivons l’équation obtenue à partir des discrétisations retenues :

22

2

x

)t,xx(u)t,x(u2)t,xx(u

x

u

ΔΔ++−Δ−=

∂∂ ;

t

)t,x(u)tt,x(u

t

u

Δ−Δ+=

∂∂  

cela donne2x

)t,xx(u)t,x(u2)t,xx(ua

t

)t,x(u)tt,x(u

ΔΔ++−Δ−

−Δ+ 

soit en mettant les termes inconnus à gauche et les termes connus à droite

( ))t,xx(u)t,x(u2)t,xx(ux

ta)t,x(u)tt,x(u

2Δ++−Δ−

ΔΔ

+=Δ+  

en posant2x

tar 

ΔΔ

=  

( ))t,xx(u)t,x(u2)t,xx(ur )t,x(u)tt,x(u Δ++−Δ−+=Δ+  

On peut aussi écrire cette équation sous la forme

0)t,xx(ru)t,xx(ru)t,x(u)r 21()tt,x(u =Δ+−Δ−−+−+Δ+  

Calqué sur la discrétisation du domaine, on obtient l’image de l’équation

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt

-1+2r 

1

-r -r 

Schéma explicite

Ce schéma est extrêmement simple à mettre en œuvre. Mais on montre qu’il est

conditionnellement stable, à savoir il faut respecter 2

1

x

tasoit

2

1r 

2≤

ΔΔ

≤ . On a donc

une contrainte sur le choix des pas de temps et d’espace qui se traduit généralement par des

 pas de temps très petits une fois le pas d'espace choisi.

On montre que pour r =1/6, l'erreur de troncature est minimalisée.

5.4.4. Le schéma impliciteEcrivons l’équation obtenue à partir des discrétisations retenues :

22

2

x

)tt,xx(u)tt,x(u2)tt,xx(u

x

u

ΔΔ+Δ++Δ+−Δ+Δ−

=∂∂

;t

)t,x(u)tt,x(u

t

u

Δ−Δ+

=∂∂

 

cela donne 2x

)tt,xx(u)tt,x(u2)tt,xx(u

at

)t,x(u)tt,x(u

Δ

Δ+Δ++Δ+−Δ+Δ−

−Δ+ 

Page 37: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 37/85

Modélisation

FN CRES

37

( ))tt,xx(u)tt,x(u2)tt,xx(ux

ta)t,x(u)tt,x(u

2Δ+Δ++Δ+−Δ+Δ−

ΔΔ

=−Δ+  

en posant2x

tar 

ΔΔ

= et en mettant les termes inconnus à gauche et les termes connus à droite

)t,x(u)tt,xx(ru)tt,x(u)r 21()tt,xx(ru =Δ+Δ+−Δ+++Δ+Δ−−  

On peut aussi écrire cette équation sous la forme

0)t,x(u)tt,xx(ru)tt,x(u)r 21()tt,xx(ru =−Δ+Δ+−Δ+++Δ+Δ−−  

Calqué sur la discrétisation du domaine, on obtient l’image de l’équation

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt 1+2r 

-1

-r -r 

Schéma implicite

Il faut donc écrire cette équation en chaque nœud inconnu. On a à résoudre ensuite ce système

d’équations.

Supposons que la dimension d’espace entre 0 et L soit divisée en N nœuds. Par ailleurs, les

nœuds marqués d’un carré noir sont connus :

-  la condition initiale à t = 0-  les conditions limites à x = 0 et x = L (conditions de Dirichlet)

0

x

t

0

2

1 2 3 N N-1 N-2i-1 i i+1

1

u = g1 u = g2

On écrit donc les équations sur une ligne quelconque au temps t en fonction de la ligne du

dessous au temps t - Δt et en mettant les termes inconnus à gauche et les termes connus à

droite.

On notera le nœud à l’abscisse i et au temps t : t

iu

nœud 0 : ce nœud est connu, on n’écrit donc pas d’équation (il est intégré dans

l’équation suivante)

nœud 1 : t

1

1t

2

1t

1

1t

0 uruu)r 21(ru =−++− +++ or  1

1t

0 gu =+ d’où

1t1

1t2

1t1 rguruu)r 21( +=−+ ++  

Page 38: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 38/85

Modélisation

FN CRES

38

nœud 2 : t

2

1t

3

1t

2

1t

1 uruu)r 21(ru =−++− +++  

nœud 3 : t

3

1t

4

1t

3

1t

2 uruu)r 21(ru =−++− +++  

…..

nœud i : t

i

1t

1i

1t

i

1t

1i uruu)r 21(ru =−++− ++

++−  

…..nœud N-2 : t

2 N

1t

1 N

1t

2 N

1t

3 N uruu)r 21(ru −+−

+−

+− =−++−  

nœud N-1 : t

1 N

1t

 N

1t

1 N

1t

2 N uruu)r 21(ru −++

−+− =−++− or  2

1t

 N gu =+ d’où

2

t

1 N

1t

 N

1t

1 N rguruu)r 21( +=−+ −++

−  

nœud N : ce nœud est connu, on n’écrit donc pas d’équation (il est intégré dans

l’équation précédente)

On peut mettre ces équations sous forme matricielle.

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

+

+

=

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

+−

−+−

−+−

−+−

−+−−+

−+−

+−

+

+

+

+

+−

+−

++++

2

t

1 N

t

2 N

t

i

t

3

t

2

1t1

1t

1 N

1t

2 N

1t

i

1t

3

1t

2

1t1

1t

1 N

1t

2 N

1t

i

1t

3

1t

2

1t

1

rgu

u

...

u

...

u

u

rgu

u

u

...

u

...

u

u

u

r 21r 

r r 21r 

.........

r r 21r 

.........

r r 21r 

r r 21r r r 21

uuuuuu

 

 Nous avons un système tri diagonal qui peut se résoudre par l’algorithme du double balayage

(voir le dernier chapitre).

On montre que le schéma implicite est inconditionnellement stable, quelque soit la valeur de

r, donc des pas d’espace et de temps.

5.4.5. Le schéma explicite à 2 pasOn n’insistera pas sur ce schéma car on montre qu’il est toujours instable.

5.4.6. Le schéma pondéréIl est intéressant d’écrire un schéma de type particulier qu’on peut qualifier de schéma mixte

explicite – implicite ou schéma pondéré.

 Nous reprendrons la notation précédente de la forme t

iu

La dérivée en temps s'écritt

uu

t

u t

i

1t

i

Δ−

=∂∂ +

 

Ecrivons la dérivée seconde par rapport à x aux 2 instants t et t+1.

Page 39: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 39/85

Modélisation

FN CRES

39

2

t

1i

t

i

t

1i

t

i

2

2

x

uu2u

x

u

Δ+−

=⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂ +− ;

2

1t

1i

1t

i

1t

1i

1t

i

2

2

x

uu2u

x

u

Δ+−

=⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂ +

+++

+

 

Et la dérivée qu'on va utiliser à une abscisse i quelconque comme une moyenne pondérée par 

un coefficient α des 2 expressions précédentes :

t

i

2

21t

i

2

2

i

2

2

x

u)1(

x

u

x

u⎟⎟

 ⎠ ⎞

⎜⎜⎝ ⎛ 

∂∂α−+⎟⎟

 ⎠ ⎞

⎜⎜⎝ ⎛ 

∂∂α=⎟⎟

 ⎠ ⎞

⎜⎜⎝ ⎛ 

∂∂ +

avec 0 ≤ α ≤ 1

L'EDP que l'on cherche à résoudre se discrétise comme ci-dessous :

⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

Δ+−

α−+Δ

+−α=

Δ− +−

++

++−

+

2

t

1i

t

i

t

1i

2

1t

1i

1t

i

1t

1i

t

i

1t

i

x

uu2u)1(

x

uu2ua

t

uu 

En posant2x

tar 

Δ

Δ=  

0uu2u)1(r uu2ur uu t

1i

t

i

t

1i

1t

1i

1t

i

1t

1i

t

i

1t

i =+−α−−+−α−− +−++

++−

+  

On obtient un schéma à 6 points

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt1+ 2r α 

-1+2r(1-α)

-r α -r α 

-r(1-α)-r(1-α)

Schéma pondéré

L'étude de la stabilité de ce schéma pondéré montre que :

•  Pour 0 ≤ α <2

1, ce schéma est stable à condition d'avoir 

)21(2

1r 

α−≤  

•  Pour 2

1 ≤ α ≤ 1, ce schéma est toujours stable

On peut retrouver les cas précédents selon la valeur de α :

•  α = 0 : on retrouve le schéma explicite précédent.

•  α = 1 : on retrouve le schéma implicite précédent; on parle de schéma implicite pur 

car il n'y a qu'un seul nœud au temps t + 1.

Enfin, la valeur de2

1=α est particulière; elle conduit au schéma appelé Crank-Nicholson qui

est stable sans condition, et dont l'erreur de troncature est la plus faible possible parmi toutes

les valeurs de α possibles.

Page 40: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 40/85

Modélisation

FN CRES

40

 

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt1+ r 

-1 + r 

-r / 2 -r / 2

-r / 2-r / 2

Schéma de Crank-Nicholson

Ce schéma est très utilisé.

5.4.7. Cas d'un schéma explicite stable

Ce schéma est présenté ici afin de montrer au lecteur que les possibilités sont nombreuses endifférences finies pour discrétiser une EDP.

Il s'agit du schéma de Durfort-Frankel.

Les différents termes de l'équation de la diffusion sont ainsi discrétisés :

2

t

1i

1t

i

1t

i

t

1i

2

2

x

uuuu

x

u

Δ+−−

=∂∂ +

−+−  

t2

uu

t

u 1t

i

1t

i

Δ−

=∂∂ −+

 

Ce qui donne :

2

t

1i

1t

i

1t

i

t

1i

1t

i

1t

i

x

uuuua

t2

uu

Δ+−−

=Δ− +

−+−

−+

 

Ce schéma est toujours stable mais il n'est consistant que si Δt est très inférieur à Δx (on a

l'erreur de troncature2

2

x

ta)u(R 

Δ

Δ≈Δ ). On parle alors de consistance conditionnelle. Ceci se

  produit généralement quand on a une discrétisation dite croisée, c'est-à-dire que la

discrétisation par rapport à une variable fait intervenir plusieurs niveaux d'une autre variable.

5.4.8. Cas d'un problème parabolique 2D

On va reprendre l'équation de la diffusion, mais à 2 dimensions d'espace ⎟⎟ ⎠

 ⎞

⎜⎜⎝ 

⎛ 

+∂

=∂

∂2

2

2

2

y

u

x

u

at

u

 

 Nous n'allons pas reprendre tous les calculs ici, mais simplement indiquer les discrétisations

utilisées et à quelle forme d'équation on aboutit.

Schéma explicite

t

uu

t

u t

i

1t

i

Δ−

=∂∂ +

;2

t

 j,1i

t

 j,i

t

 j,1i

2

2

x

uu2u

x

u

Δ

+−=

∂∂ +−

;2

t

1 j,i

t

 j,i

t

1 j,i

2

2

y

uu2u

y

u

Δ

+−=

∂∂ +−

 

Schéma implicite pur 

tuu

tu

t

i

1t

i

Δ−=

∂∂

+

;2

1t

 j,1i

1t

 j,i

1t

 j,1i

2

2

xuu2u

xu

Δ+−=

∂∂

+

+

++

− ;2

1t

1 j,i

1t

 j,i

1t

1 j,i

2

2

yuu2u

yu

Δ+−=

∂∂

+

+

++

−  

Page 41: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 41/85

Modélisation

FN CRES

41

 

Schéma pondéré

t

uu

t

u t

i

1t

i

Δ−

=∂∂ +

;

2

t

 j,1i

t

 j,i

t

 j,1i2

1t

 j,1i

1t

 j,i

1t

 j,1i2

2

x

uu2u)1(

x

uu2u

xu Δ

+−α−+Δ

+−α=∂∂

+−

+

+

++

− avec 0 ≤ α ≤ 1

2

t

1 j,i

t

 j,i

t

1 j,i

2

1t

1 j,i

1t

 j,i

1t

1 j,i

2

2

y

uu2u)1(

y

uu2u

y

u

Δ

+−α−+

Δ

+−α=

∂∂ +−

++

++−

avec 0 ≤ α ≤ 1

Les nœuds concernés par chaque schéma sont indiqués ci-dessous.

xt

t+Δt

y

Schéma

explicite pur 

Schéma

implicite pur Schéma

Pondéré

α 

(1 - α)

En posant2x

x

tar 

ΔΔ

= et2y

y

tar 

ΔΔ

=  

•  Pour 0 ≤ α <2

1, ce schéma est stable à condition d'avoir 

)21(2

1r r  yx α−

≤+  

•  Pour 2

1 ≤ α ≤ 1, ce schéma est toujours stable

Par ailleurs :

•  α = 0 : on retrouve le schéma explicite pur.

•  α = 1 : on retrouve le schéma implicite pur.

•  α =2

1: c'est le schéma de Crank-Nicholson.

5.5. Cas d'une EDP hyperbolique (évolutive)

 Nous allons prendre comme exemple l'équation dite de convection 0x

uc

t

u=

∂∂

+∂∂

.

Elle correspond par exemple à un transport de contaminant en rivière sans diffusion où c est la

vitesse de l'eau et peut s'établir facilement en établissant un bilan de polluant comme au début

de ce cours.

Page 42: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 42/85

Modélisation

FN CRES

42

5.5.1. Schémas explicites

La dérivée par rapport au temps est toujours discrétisée selon le schémat

uu

t

u t

i

1t

i

Δ−

=∂∂ +

.

Il y a ensuite différents cas selon la façon de discrétiser la dérivée en espace.

•  xuu

xu

t

1i

t

i

Δ−=∂∂−  

Ce qui conduit à l'équation discrétisée 0x

uuc

t

uu t

1i

t

i

t

i

1t

i =Δ−

− −+

 

En posantx

tcr 

ΔΔ

= , nous avons ( ) 0uur uu t

1i

t

i

t

i

1t

i =−+− −+ .

Le schéma obtenu est nommé explicite décentré amont

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt

1

r - 1-r 

Schéma explicite décentré amont 

Ce schéma est stable si c > 0 et r ≤ 1, instable si c < 0

• x

uu

x

u t

i

t

1i

Δ−

=∂∂ +  

Ce qui conduit à l'équation discrétisée 0x

uuc

t

uu t

i

t

1i

t

i

1t

i =Δ

−+

Δ− +

+

 

En posantx

tcr 

ΔΔ

= , nous avons ( ) 0uur uu t

i

t

1i

t

i

1t

i =−+− ++ .

Le schéma obtenu est nommé explicite décentré aval

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt

1

- r - 1 r 

Schéma explicite décentré aval

Ce schéma est stable si c < 0 et r  ≤ 1, instable si C > 0

Page 43: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 43/85

Modélisation

FN CRES

43

• x2

uu

x

u t

1i

t

1i

Δ−

=∂∂ −+  

Ce qui conduit à l'équation discrétisée 0x2

uuc

t

uu t

1i

t

1i

t

i

1t

i =Δ−

− +++

 

En posantxtcr 

ΔΔ= , nous avons ( ) 0uur uu t1i

t1i

ti

1ti =−+− −++ .

Le schéma obtenu est nommé explicite centré

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt

1

-1 r/2-r/2

Schéma explicite centré

Ce schéma est toujours instable

5.5.2. Schémas implicites

La dérivée par rapport au temps est toujours discrétisée selon le schéma

t

uu

t

u t

i

1t

i

Δ

−=

∂ +

.

Il y a ensuite différents cas selon la façon de discrétiser la dérivée en espace.

• x

uu

x

u 1t

1i

1t

i

Δ−

=∂∂ +

−+

 

Ce qui conduit à l'équation discrétisée 0x

uuc

t

uu 1t

1i

1t

i

t

i

1t

i =Δ−

− +−

++

 

En posantx

tcr 

ΔΔ

= , nous avons ( ) 0uur uu 1t

1i

1t

i

t

i

1t

i =−+− +−

++ .

Le schéma obtenu est nommé implicite décentré amont

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt

-1

r + 1-r 

Schéma implicite décentré amontCe schéma est stable si c > 0 , instable si c < 0

Page 44: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 44/85

Modélisation

FN CRES

44

 

• x

uu

x

u 1t

i

1t

1i

Δ−

=∂∂ ++

+  

Ce qui conduit à l'équation discrétisée 0

x

uuc

t

uu 1t

i

1t

1i

t

i

1t

i =

Δ

−+

Δ

− +++

+

 

En posantx

tcr 

ΔΔ

= , nous avons ( ) 0uur uu 1t

i

1t

1i

t

i

1t

i =−+− +++

+ .

Le schéma obtenu est nommé implicite décentré aval

0 L

xt0 

t

x x+Δxx-Δx

t

t+Δt1 - r 

- 1

Schéma implicite décentré aval

Ce schéma est stable si c < 0, instable si c > 0

• x2

uu

x

u 1t

1i

1t

1i

Δ−

=∂∂ +

−++  

Ce qui conduit à l'équation discrétisée 0x2uuc

tuu

1t

1i

1t

1i

t

i

1t

i =Δ−+Δ−+

+

+

+

+

 

En posantx

tcr 

ΔΔ

= , nous avons ( ) 0uu2

r uu 1t

1i

1t

1i

t

i

1t

i =−+− +−

++

+ .

Le schéma obtenu est nommé implicite centré

0 L

xt0 

t

x x+Δxx-Δx

tt+Δ

t

1

-1

2r −   2r   

Schéma implicite centré

Ce schéma est stable si r  ≤ 1

Page 45: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 45/85

Modélisation

FN CRES

45

5.6. Maillages particuliersOn indique dans ce paragraphe comment opérer les discrétisation quand on a des pas d'espace

irréguliers. Les pas de temps irréguliers ne posent pas de problème particulier car le pas de

temps peut être choisi à chaque itération pour passer d'un pas de temps à l'autre.

5.6.1. Maillage irrégulier en xOn traite le cas de l'axe Ox, mais on retrouve le même principe quelque soit l'axe. On se met

donc dans le cas suivant:

u0 u2u1

h1 h2u(x-h1) u(x) u(x+h2) x

Et on cherche une discrétisation des dérivées première et seconde. Nous avons :

...x

u

!4

h

x

u

!3

h

x

u

!2

h

x

uh)x(u)hx(uu

4

44

2

3

33

2

2

22

2222 +

∂∂

+∂∂

+∂∂

+∂∂

+=+=  

...x

u

!4

h

x

u

!3

h

x

u

!2

h

x

uh)x(u)hx(uu

4

441

3

331

2

221

111 +∂∂+

∂∂−

∂∂+

∂∂−=−=  

5.6.1.1   Discrétisation de la dérivée première

D’après la première équation

)h(oh

)x(u)hx(u

x

u2

2

2 +−+

=∂∂

: discrétisation avant

D’après la deuxième équation

)h(oh

)hx(u)x(u

x

u1

1

1 +−−

=∂∂

: discrétisation arrière

En multipliant la première équation par  2

1h , la deuxième par  2

2h et en soustrayant :

...x

u)hh(hhu)hh)(hh(...

x

u)hhhh(u)hh(uhuh 2121021211

2

22

2

10

2

2

2

11

2

22

2

1 +∂∂

++−+=+∂∂

++−=−

...u)hh(hh

hu

)hh(hh

hu

)hh(hh

)hh)(hh(

x

u1

2121

2

22

2121

2

10

2121

2121 ++

−+

++

−+−=

∂∂

 

...uh

hu

h

h

)hh(

1u

h

1

h

1

x

u1

1

22

2

1

21

0

21

+⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ −

++⎟⎟

 ⎠

 ⎞⎜⎜⎝ 

⎛ −=

∂∂

: discrétisation centrée

 Remarque : si 21 hh = = h, alors ...h2

uu

x

u 12 +−

=∂∂

 

5.6.1.2   Discrétisation de la dérivée seconde

En multipliant la première expression par  1h , la deuxième par  2h et en additionnant :

Page 46: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 46/85

Modélisation

FN CRES

46

...x

u)hhhh(

2

1u)hh(uhuh

2

2

2

2

11

2

20211221 +∂

∂+++=−  

...u)hh(hh

hh2

)hh(hh

uhuh2

x

u0

2121

21

2121

1221

2

2

++

+−

++

=∂∂

 

...uhh

2

h

u

h

u

hh

2

x

u0

212

2

1

1

21

2

2

+−⎟⎟ ⎠ ⎞

⎜⎜⎝ ⎛  +

+=

∂∂ : discrétisation centrée

 Remarque : si 21 hh = = h, alors  ...h

uu2u...u

h

2

h

u

h

u

h2

2

x

u2

20102

21

2

2

++−

=+−⎟ ⎠

 ⎞⎜⎝ 

⎛  +=∂∂

qui est la

discrétisation qu'on a déjà utilisée abondamment.

5.6.2. Maillage différent en x et en y

On reprend le premier exemple traité (équation de Laplace dans un domaine carré). On enreste donc à un problème 2D, l'extension à un problème 3D ne posant pas de problème

 particulier.

On va simplement étudier la discrétisation en considérant Δx ≠ Δy

Les discrétisations des 2 dérivées secondes donnent :

)x(ox

uu2u

x

u 2

2

 j,1i j,i j,1i

2

2

Δ+Δ

+−=

∂ +− 

)y(oy

uu2u

y

u 2

2

1 j,i j,i1 j,i

2

2

Δ+Δ

+−=

∂ +− 

En sommant ces 2 expressions

)y(o)x(oy

uu2u

x

uu2u

y

u

x

u 22

2

1 j,i j,i1 j,i

2

 j,1i j,i j,1i

2

2

2

2

Δ+Δ+Δ

+−+

Δ

+−=

∂+

∂ +−+− 

)y(o)x(oy

u

y

u

y

2

x

2u

x

u

x

u

y

u

x

u 22

2

1 j,i

2

1 j,i

22 j,i2

 j,1i

2

 j,1i

2

2

2

2

Δ+Δ+Δ

+⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

Δ+

Δ−

Δ+

Δ=

∂∂

+∂∂ +−+−

 

En calquant cette expression sur le maillage :

Ui,j

Ui,j-1

Ui,j+1

Ui+1,jUi-1,j 22 y

2

x

2

Δ−

Δ−  

2y1

Δ 

2y1

Δ 

2x1

Δ  2x

 

On aboutit donc à des expressions un peu plus compliquée, mais ne présentant pas de

difficultés particulières.

Page 47: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 47/85

Modélisation

FN CRES

47

 

5.7. Traitement des termes non linéaires

Très souvent, les coefficients qui interviennent dans les EDP sont eux-mêmes fonctions de la

fonction cherchée.

Par exemple, dans l'équation de la diffusion 1D, on a2

2

x

u)u(a

t

u

∂∂=

∂∂ .

La discrétisation de cette équation pose alors un problème dès qu'on utilise un schéma autre

que le schéma explicite (ce qui est souvent le cas). En effet, on travaille entre les pas de temps

(t) et (t+1), et la valeur de a(u) qu'il faut considérer est alors une valeur entre )u(a t et

)u(a 1t+ , or on ne connait pas 1tu + , car c'est ce que l'on cherche.

Pour passer du temps (t) au temps (t+1), on va donc procéder à une série d'itérations visant à

obtenir une bonne valeur de a(u) sur l'intervalle de temps étudié. Il faut bien noter que ces

itérations se superposent au système d'itérations précédent.

t + 2t + 1t

 Itérations pour passer d'un

 pas de temps à l'autre

 Itérations supplémentaires à

l'intérieur d'un pas de temps

temps

Considérons le schéma le plus général qu'on ait écrit, le schéma pondéré, où la dérivée

seconde est discrétisée par 

t

i

2

21t

i

2

2

i

2

2

x

u)1(

x

u

x

u⎟⎟

 ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂

α−+⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂

α=⎟⎟ ⎠

 ⎞⎜⎜⎝ 

⎛ 

∂∂

+

avec 0 ≤ α ≤ 1

Pour la première itération (k = 1), on commence par évaluer une valeur de u au pas de temps

t+1 (notée 1t

1u + ). Ceci peut être réalisé en supposant une variation linéaire de u :1tt1t

1 uu2u −+ −=  

Ce qui permet de calculer une valeur de a1(u), dont on peut prendre la moyenne pondérée sur 

les 2 pas de temps

)u(a)1()u(a)u(a t1t

11 α−+α= +  

Ceci permet de calculer une valeur de 1t

2u + par résolution du système matriciel.

Et on recommence les itérations

)u(a)1()u(a)u(a t1t

k k  α−+α= +  

  jusqu'à convergence de la solution obtenue pour ut+1, c'est-à-dire quand la valeur de ut+1 ne

 bouge plus avec une certaine précision fixée.

Page 48: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 48/85

Modélisation

FN CRES

48

 

6. LA METHODE DES ELEMENTS FINIS

6.1. Comparaison Différences finies – Eléments finisLe schéma ci-dessous compare les 2 méthodes en termes de démarche. La méthode des

différences finies cherche à résoudre exactement une équation approchée, alors que laméthode des éléments finis cherche une solution approchée à une équation exacte.

Equation physique continue

aux Dérivées Partielles

reformulation discrète approchée reformulation intégrale exacte

de l’équation d’origine de l’équation d’origine

résolution exacte de résolution approchée de

l’équation approchée l’équation exacte

Méthodes d’approximation Méthodes d’approximationd’équation de solution

(Différences finis) (Eléments finis)

6.2. Equivalence de problèmes

6.2.1. ObjetL’objet de ce paragraphe est de montrer comment on peut avoir des problèmes d’expression

totalement différente mais dont la solution est identique. En effet, la méthode des éléments

finis suppose qu’on transforme l’EDP en un autre problème ayant même solution; et c’est ce

nouveau problème qu’on résoudra.

Bien sûr, la mise en place de ces équivalences de problèmes va impliquer certaines conditions

sur la formulation des problèmes eux-mêmes.

6.2.2. Equivalence entre système matriciel et problème de minimisationSoit un système d’équations linéaires (n équations, n inconnues) de la forme

Au = b avec A = [aij] matrice carrée de dimension (n,n);

u = [ui] vecteur de dimension (n), solution du système;

 b = [bi] vecteur de dimension (n);

On montre que si la matrice A est symétrique, la solution u du système matriciel est aussi la

solution d’un problème de minimum qui s’exprime ainsi :

trouver u telle que la quantité I(u) soit minimum avec

Page 49: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 49/85

Modélisation

FN CRES

49

I(u) = [ ]u u

a a

a

a a

u

u

n

n

ij

n nn n

1

11 1

1

1

..

..

.. ..

..

..

⎢⎢⎢

⎥⎥⎥

⎢⎢⎢

⎥⎥⎥

- 2 [ ] b b

u

u

n

n

1

1

.. ..

⎢⎢⎢

⎥⎥⎥

= (ut A u ) - 2 (bt u)

qu’on peut noter I(u) = a(u,u) - 2 L(u) où a(u,u) = u

t

A u

 

L(u) = bt u

remarque : la quantité (ut A u) est appelée ‘forme quadratique’ associée à la matrice A.

En effet, I(u) = a u u b uij i j j

n

i

n

i ii

n

== =∑∑ ∑−

11 1

2

La condition pour que u soit le vecteur minimisant I(u) s’écrit : ∀ p=1,n :∂∂I u

u p

( )= 0

∂∂I u

ua u a u b

 p pj j

 j

n

ip ii

n

 p

( )= + −

= =∑ ∑

1 1

2 = a u a u b pi ii

n

ip ii

n

 p= =∑ ∑+ −

1 1

2 0 ∀ p=1,n

On utilise ici le fait que la matrice A est symétrique : aip = a pi d’où

∂∂I u

ua u b

 p

 pi i

i

n

 p

( )= −

=∑2 2

1

= 0 → a u b pi i

i

n

 p=

∑ =1

  ∀  p=1,n soit Au = b -cqfd -

Ainsi, ce problème matriciel et ce problème de minimum ont la même solution.

6.2.3. Equivalence avec un problème de minimumSoit le problème de minimum précédent noté sous le forme I(u) = a(u,u) - 2L(u).

On montre que si a est une forme bilinéaire, symétrique et positive et L une forme linéaire

{   forme linéaire : soit E un espace vectoriel, une forme linéaire dans E est uneapplication linéaire de E dans R.

 forme bilinéaire :linéaire en u : a(λ1u1 + λ2u2 , v) = λ1a(u1,v) + λ2a(u2,v)

linéaire en v : a(u, μ1v1 + μ1v2) = μ1a(u,v1) + μ2a(u,v2)

symétrique : a(u,v) = a(v,u)

 positive : a(u,u) ≥ 0 }

alors le problème de minimum est équivalent à :

trouver u∈V tel que ∀ v∈V : a(u,v) = L(v).

V peut être défini comme l’espace vectoriel où on exprime la solution u. Plus loin, on verra

apparaître des conditions sur cet espace.En effet : en se plaçant dans V

Si u vérifie I(u) minimum, alors on doit avoir : ∀ v∈V, ∀ ε∈R : I(u + εv) ≥ I(u)

soit I(u + εv) = a(u + εv , u + εv) - 2 L(u + εv)

en utilisant la bilinéarité de a et la linéarité de L :

I(u + εv) = a(u,u) + εa(u,v) + εa(v,u) + ε²a(v,v) - 2 L(u) - 2 εL(v)

en utilisant la symétrie de a :

I(u + εv) = I(u) + 2ε [a(u,v) - L(v)] + ε²a(v,v)

or si on veut I(u + εv) ≥ I(u) cela implique 2ε [a(u,v) - L(v)] + ε²a(v,v) ≥ 0

en utilisant le fait que a est positive : a(v,v) ≥ 0, pour que l’expression ci-dessus soit

toujours positive ∀ ε, il faut avoir : a(u,v) - L(v) = 0 soit a(u,v) = L(v) -cqfd -

Ainsi, avec les conditions émises sur a et L , les deux problèmes ont la même solution.

Page 50: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 50/85

Modélisation

FN CRES

50

L’équation a(u,v) = L(v) est appelée équation d’Euler associée au problème de minimum.

6.2.4. Récapitulatif Le schéma suivant récapitule les équivalences de problèmes précédentes

Trouver u telle que

Au = b

Trouver u telle que

I(u) = a(u,u) – 2L(u)

soit minimum

Système matriciel

Problème de minimum

Si A symétrique

Si a(u,v) forme bilinéaire, symétrique, positive

 L(u) forme linéaire

Equation d'Euler associée Soit V en espace vectoriel où on exprime la solution

Trouver v ∈ V tel que

∀ u ∈ V a(u,v) = L(v)

Ainsi, on a une série d’équivalence de problèmes qui s’enchainent ainsi :

système matriciel → problème de minimum → équation d’Euler associée

et ces équivalences sont liées à certaines conditions.

Ceci est à la base de la méthode dite de Ritz, la méthode des éléments finis consistant

finalement en un choix très judicieux d'un sous-espace vectoriel de V.

6.2.5. Application aux EDPOn va voir ici le principe de transformation d’une EDP en un problème de minimum.

Soit G un opérateur différentiel et l’équation en u : Gu = f sur un domaine Ω de frontière Γ.

On ne change pas l’équation si on multiplie à gauche et à droite par un vecteur v∈V, soit

v Gu = vf 

On peut intégrer cette équation (une condition sur v, donc sur V, est qu’il soit intégrable) :

( ) ( )v Gu d vf dΩ ΩΩ Ω∫ ∫ =  

Si on peut maintenant ramener l’expression de ( )v Gu dΩΩ∫  à une forme bilinéaire,

symétrique et positive et ( )vf dΩΩ∫  à une forme linéaire, alors le problème reviendra à

a(u,v) = L(v)

soit un problème équivalent au problème de minimum précédent.

Cette dernière expression est appelée ‘forme variationnelle’ de l’EDP.

(On verra plus loin concrètement comment on passe de l’EDP à sa forme variationnelle)

 Remarque 1:

 L’espace V où on exprime la solution est l’espace L² des fonctions dites de carré sommable.Ces fonctions f (qui peuvent être complexes) sont telles que | ( )|²f x dx

−∞

+∞

∫  existe. Le produit de

Page 51: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 51/85

Modélisation

FN CRES

51

2 fonctions appartenant à L² appartient à L1

, espace des fonctions sommables, c’est-à-dire

intégrables au sens de Lebesgue (fonctions f telles que f x( )dx−∞

+∞

∫  existe).

 L² et L1

sont des espaces vectoriels.

Par ailleurs, on démontre qu’une fonction bornée f, nulle en dehors d’un intervalle fini (a,b),

est sommable. Si, de plus, f est intégrable au sens de Riemann sur (a,b), alors son intégrale de Lebesgue sur cet intervalle est égale à son intégrale de Riemann.

On verra qu’effectivement l’espace vectoriel où on va exprimer la solution de l’EDP est un

espace de fonctions nulles en dehors du domaine d’étude Ω  de l’EDP et qu’il suffira donc

d’intégrer les équations au sens de Riemann (pourvu qu’on sache l’intégrer).

Ce que nous venons de voir dans R s’étend aussi à Rn , y compris donc à un domaine Ω où on

cherche la solution d’une EDP.

On parlera donc de L²(Ω ) espace des fonctions de carré sommable dans Ω .

6.2.6. Approximation interne

La solution théorique, lorsqu’elle existe, du problème a(u,v)=L(v) s’exprime généralement

dans un espace vectoriel V de dimension infinie. Elle s’écrit sous la forme

u x N xi ii

( ) ( )r r

==

∑α1

 

où N xi ( )r

forme une base de cet espace vectoriel.

Réaliser une approximation de cette solution consiste alors à exprimer la solution dans un

sous-espace vectoriel noté Vh de V de dimension finie : c’est l’approximation interne.

Choisissons une base finie de vecteurs Ni, i=1,n (on omettra, pour la suite de ce paragraphe,

les composantesrx des vecteurs) engendrant le sous-espace vectoriel Vh.

Le problème s’exprime désormais ainsi :

Trouver uh∈Vh tel que ∀vh∈Vh : a(uh, vh) = L(vh).

vh peut s’exprimer dans la base des Ni : v Nh j j j

n

==∑β

1

; en remplaçant vh :

a(uh , β j j j

n

 N=∑

1

) = L( β j j j

n

 N=∑

1

) en utilisant la linéarité de L :

β β j h j j

n

 j j j

n

a u N L N( , ) ( )= =∑ ∑=

1 1

ou encore

β j h j j

n

 ja u N L N[ ( , ) ( )]=∑ − =1

0 ceci devant être vrai ∀vh, donc ∀β j :

→  ∀  j =1,n : a(uh, N j) = L(N j) d’où le nouveau problème :

Trouver uh∈Vh tel que ∀ j =1,n : a(uh, N j) = L(N j).

uh peut aussi s’exprimer dans la base des Ni : u Nh i ii

n

==

∑α1

; en remplaçant uh :

a( α i i j

n

 N=∑

1

, N j ) = L(N j) en utilisant la linéarité de a :

α i i ji

n

 ja N N L N( , ) ( )=∑ =

1

d’où le problème :

Page 52: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 52/85

Modélisation

FN CRES

52

Trouver  αi tel que ∀  j =1,n : α i i ji

n

 ja N N L N( , ) ( )=∑ =

1

 

ce qui est un système de n équations à n inconnues que l’on peut écrire :

a N N a N N a N N

a N N a N N a N N

a N N a N N a N N

L N

L N

L N

n

n

n n n n n n

( , ) ( , ) ... ( , )

( , ) ( , ) ... ( , )

... ...

( , ) ( , ) ... ( , )

...

( )

( )

...

( )

1 1 1 2 1

2 1 2 2 2

1 2

1

2

1

2

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

=

⎢⎢⎢⎢

⎥⎥⎥⎥

α

α

α

 

Remarques : la matrice carrée est symétrique.

les αi n’ont aucune signification physique à priori

Ainsi, on a vu qu’on avait une série d’équivalence de problèmes qui s’enchainaient ainsi :

système matriciel → problème de minimum → équation d’Euler associée → EDP

et on vient de voir que l’on pouvait remonter la filière pour, à partir d’une EDP, aboutir à unsystème d’équations dont la solution est une solution approchée de l’EDP. Ceci est à la base

de la méthode dite de Ritz, la méthode des éléments finis consistant finalement en un choix

très judicieux des vecteurs Ni, base du sous-espace vectoriel d’approximation.

6.3. Construction pratique du problème variationnel

6.3.1. Cas d’une équation différentielle

On traitera dans ce paragraphe l’exemple suivant :-u"(x) + u(x) = f(x) = 2x -1 dans l’intervalle (a,b)

on prendra toujours u’(b)=0 et on fera varier la condition en (a).

6.3.1.1  Cas d’une condition de Neumann homogènesoit u’(a) = 0

En reprenant la démarche précédente, on écrit :

∀ v∈V = { v : v∈L²(a,b) } v (-u" + u) = v f 

( ( ( )dxv vf a

 b

a

 b

− =∫ ∫ u"+u) )dx

( ( )dx ( )dx− + =∫ ∫ ∫ vu" vu vf  a

 b

a

 b

a

 b

)dx

en intégrant le premier terme par parties (possible si v’∈L²(a,b)) :

-v(b)u’(b) + v(a)u’(a) + ( ' ( )dx ( )dxv vu vf  a

 b

a

 b

a

 b

u')dx∫ ∫ ∫ + =  

comme u’(a) = u’(b) = 0 :

( ' ( )dxv vf a

 b

a

 b

u' + vu)dx∫ ∫ =  

Le premier terme est bien une forme bilinéaire, symétrique et positive, le deuxième

une forme linéaire, d’où le nouveau problème :

Page 53: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 53/85

Modélisation

FN CRES

53

trouver u∈V tel que ∀v∈V, on ait a(u,v) = L(v)

avec V = { v : v∈L²(a,b); v’∈L²(a,b) }

a(u,v) = ( 'va

 b

u' + vu)dx∫  et L(v) = ( )dxvf a

 b

∫   

6.3.1.2  Cas d’une condition de Dirichlet homogènesoit u(a) = 0

Les calculs sont identiques jusqu’à :

-v(b)u’(b) + v(a)u’(a) + ( ' ( )dx ( )dxv vu vf  a

 b

a

 b

a

 b

u')dx∫ ∫ ∫ + =  

Ici, u’(b) = 0 d’où

v(a)u’(a) + ( ' ( )dxv vf a

 b

a

 b

u' + vu)dx∫ ∫ =  

Pour ramener le premier terme à une forme bilinéaire, symétrique et positive, on va se servir 

de la liberté que l’on a de choisir l’espace vectoriel où on va exprimer la solution. Ici, on peut

le choisit tel que v(a) = 0.D’où le nouveau problème :

trouver u∈V tel que ∀v∈V, on ait a(u,v) = L(v)

avec V = { v : v∈L²(a,b); v’∈L²(a,b); v(a)=0 }

a(u,v) = ( 'va

 b

u' + vu)dx∫  et L(v) = ( )dxvf a

 b

∫   

6.3.1.3  Cas d’une condition de Neumann non homogènesoit u’(a) = g N 

Les calculs sont identiques jusqu’à :-v(b)u’(b) + v(a)u’(a) + ( ' ( )dx ( )dxv vu vf  

a

 b

a

 b

a

 b

u')dx∫ ∫ ∫ + =  

Comme u’(a)=g N et u’(b)=0, nous avons :

( ' ( )dxv vf a

 b

a

 b

u' + vu)dx∫ ∫ = - v(a)g N

Le vecteur v étant choisit (par le choix de la base où on exprime la solution - voir 6-2.5.) et g N 

étant connu, il n’y a pas de difficulté à calculer le terme v(a)g N 

D’où le nouveau problème :

trouver u∈V tel que ∀v∈V, on ait a(u,v) = L(v)

avec V = { v : v∈L²(a,b); v’∈L²(a,b)}a(u,v) = ( 'v

a

 b

u' + vu)dx∫  et L(v) = ( )dxvf a

 b

∫  - v(a)g N

6.3.1.4  Cas d’une condition de Dirichlet non homogènesoit u(a) = gD 

Les calculs sont identiques jusqu’à :

-v(b)u’(b) + v(a)u’(a) + ( ' ( )dx ( )dxv vu vf  a

 b

a

 b

a

 b

u')dx∫ ∫ ∫ + =  

En choisissant v tel que v(a)=0, nous avons :

( ' ( )dxv vf a

 b

a

 b

u' + vu)dx∫ ∫ =  

Page 54: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 54/85

Modélisation

FN CRES

54

avec V = { v : v∈L²(a,b); v’∈L²(a,b); v(a)=0 } 

Le problème qui apparaît alors est que u∉V puisque u(a)=gD et il faudrait u(a)=0.

On définit un nouvel espace vectoriel U = { u : u∈L²(a,b); u’∈L²(a,b); u(a)=gD }

Tout vecteur de U peut s’écrire : u = u0 [∈U] + w [∈V]

  puisque u0(a)=gD et w(a)=0, on retrouve bien u(a)=gD 

Le nouveau problème peut alors s’écrire :

Soit u = u0 + w la solution de l’équation ; u0∈U0 

trouver w∈V tel que ∀v∈V, on ait a(w,v) = L(v)

avec U0 = { u0 : u0∈L²(a,b); u0’∈L²(a,b); u0(a)=gD }

V = { v : v∈L²(a,b); v’∈L²(a,b); v(a)=0 }

a(w,v) = ( 'va

 bw' + vw)dx∫  et L(v) = ( )dxvf 

a

 b

∫   

6.3.2. Cas d’une EDPOn traitera là aussi d’un exemple qui est l’équivalent du problème précédent, mais à plusieurs

dimensions :

− + =Δu x u x f x( ) ( ) ( )r r r

sur le domaine Ω ayant une frontière Γ;

avec la condition aux limites u x( ) /r

Γ = 0 (Dirichlet homogène);

( )rx est un vecteur de composantes (x1, x2, ..., xn);

On suppose par ailleurs que Γ est suffisamment régulière pour que la normale

existe partout.

Le processus est identique au cas précédent : on multiplie par un vecteur (v) à droite et à

gauche de l’équation, et on intègre sur Ω (en notant u u x et v v x= =( ) ( )r r ) :

(( )v )d ( )d

(( )v )d ( )d ( )d

− + =

− + =

∫ ∫ ∫ ∫ ∫ 

Δ Ω Ω

Δ Ω Ω Ω

Ω Ω

Ω Ω Ω

u u fv

u uv fv 

Il faut transformer le premier terme du premier membre de l’égalité afin d’aboutir sur ce

 premier membre à une forme bilinéaire, symétrique et positive. Ceci se fait grâce à la formule

de Green, dont l’application à un problème à une dimension n’est rien d’autre que la fameuse

intégration par parties. La formule de Green s’écrit :

uv

xd v

u

xd uvn d

i ixi

∂Ω Ω Γ

Ω Ω Γ∫ ∫ ∫ + =  

où xi est la ième composante derx et nxi la ième composante de

rn , normale à Γ.

En remplaçant (u) par sa dérivée dans la formule de Green :

∂∂

∂∂

∂∂

u

x

v

xd v

u

xd v

u

xn d

i i i i

xiΩ Ω ΓΩ Ω Γ∫ ∫ ∫ + =

2

Cette relation est vraie ∀i; en sommant ces équations :

( . ) ( . )d ( ( ))d ( ( . ))dr r r r∇ ∇ + = = ∇∫ ∫  ∑∫ ∫ 

=

u v d v u v nu

xv n uxi

ii

n

Ω Δ Ω Γ ΓΩ Ω Γ Γ

∂∂1

 

Page 55: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 55/85

Modélisation

FN CRES

55

En revenant à notre exemple et en remplaçant (( ) )d−∫  Δ ΩΩ

u v obtenu par la formule de Green

( . ) ( ( . ))d ( )d ( )dr r r r∇ ∇ − ∇ + =∫ ∫ ∫ ∫  u v d v n u uv fvΩ Γ Ω Ω

Ω Γ Ω Ω 

Il suffit alors de choisir v tel que v/Γ=0 pour obtenir 

( . ) ( )d ( )dr r

∇ ∇ + =∫ ∫ ∫ u v d uv fvΩ Ω ΩΩ Ω Ω 

d’où le nouveau problème :

trouver u∈V tel que ∀v∈V, on ait a(u,v)=L(v)

avec V = { v : v∈L²(Ω);∂∂

v

x i

∈L²(Ω); v/Γ=0 }

a(u,v) = ( . )r r∇ ∇ +∫  u v uv dΩ

Ω 

L(v) = ( )dfv ΩΩ∫   

 Remarque : V est l’espace de Sobolev noté H 01  

6.4. Mise en œuvre de la méthode des éléments finis

Cette méthode est particulièrement adaptée aux problèmes dont les frontières sont d’une

géométrie irrégulière.

On découpe le domaine d’étude en éléments de géométrie simple, mais qui vont pouvoir 

couvrir la zone d’étude de façon complète (sans ‘trou’), et en suivant au mieux la frontière.

Par exemple, si le domaine d’étude Ω ⊂ R², avec une frontière Γ quelconque, on peut choisir 

des éléments Ωh triangulaires pour recouvrir Ω et en affinant autant que nécessaire au niveau

de Γ.

Ω h

Ω

 

Ce principe est généralisable dans R 3 avec des tétraèdres.

On va alors définir une formulation variationnelle sur  Ωh et une approximation interne dans

un espace vectoriel facilement manipulable sur l’ensemble des Ωh.

6.5. Cas d’une équation différentielle - Fonction linéaire parmorceaux

6.5.1. Cas de condition de Dirichlet homogèneOn reprend l'exemple précédent et les conditions limites suivantes :

-u"(x) + u(x) = f(x) = 2x - 1 sur le domaine Ω = [0,1]

avec u(0) = 0 et u(1) = 0

Page 56: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 56/85

Modélisation

FN CRES

56

 

On a vu que ce problème a la même solution que le problème suivant :

trouver u∈V tel que ∀ v∈V , on ait : a(u,v) = L(v)

avec V = { v : v∈L²(Ω) ; v’∈L²(Ω) ; v(0)=0; v(1)=0}

a(u,v) = ( ' ' )uv u v dx+∫ 01  

L(v) = ( )fv dx0

1

∫   

Choix des éléments

Le domaine d’étude est séparé en éléments de largeur (h). Cette largeur peut être différente

 pour chaque élément (dans le cas où on veut affiner la connaissance de la solution à certains

endroits)

Dans cet exemple et pour simplifier les écritures qui vont suivre, on choisit de séparer le

domaine d'étude en 5 éléments de même longueur.

x6x1 x2 x3 x4 x5 x6

0 1

ΩΩh

 

Approximation linéaire par morceaux

On va approximer la solution sur chaque élément par un morceau de droite

u 5

u4

u6

u 3 

x

u

u1  

u 2 

x 1 x 2 x 3 x 4 x 5 x6

 

 Remarque :

- il y a continuité de la fonction par morceaux

- il n'y a pas continuité de la dérivée

On va exprimer la solution u(x) cherchée dans un espace vectoriel défini par une base de

vecteurs que l’on note Ni(x). Cette base sera finie afin qu’on puisse effectivement calculer 

u(x).

u x N xi ii

n

( ) ( )==∑α

1

 

Page 57: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 57/85

Modélisation

FN CRES

57

Si la base est définie, il reste à calculer les αi .

Deux questions apparaissent à propos de la base de vecteurs Ni(x) :

- quels vecteurs ?

- combien de vecteurs ?

Il serait intéressant que les αi que l’on cherche correspondent aux ui de la fonction tels qu’ilssont définis sur le graphique précédent. Trouver les αi conduirait alors à avoir directement

une estimation de la fonction aux abscisses xi.

Comment choisir les Ni(x) pour respecter l’alinéa précédent ?

Ecrivons la solution u(x) en une abscisse particulière, par exemple x2, correspondant à une

limite d’élément et prenons comme nombre de fonctions Ni(x) le nombre d’inconnues auquel

on est confronté (6 ici)

u(x2) = α1 N1 (x2) + α2 N2 (x2) + ... + α6 N6(x2) = u2 

Pour que α2

= u2

à l’abscisse x2, il suffit de choisir les N

i(x) tels que :

 N1 (x2) = N3 (x2) = ... = N6 (x2) = 0 et N2 (x2) = 1

Pour une abscisse (x) quelconque comprise entre xi et xi+1 , u(x) est une valeur comprise entre

ui et ui+1 qui est la moyenne pondérée par les distances respectives entre (xi et x) et (xi+1 et x).

Considérons un élément Ωh et écrivons symboliquement la solution :

 Ni-1(x)= ui-1  + ui

xi-1 xi xi-1 xi xi-1 xi

u(x)

 Ni(x)ui-1 

ui  1

0

1

0

toutes les autres

fonctions N j(x)

[j≠i-1, j≠i] nullessur cet élément

+

 Afin d’identifier complètement la fonction Ni(x), considérons l’élément contigu au précédent :

= ui  + ui+1

xi xi+1 xi xi+1 xi xi+1

u(x) Ni(x) Ni+1(x

ui 

ui+1  1

0

1

0

toutes les autres

fonctions N j(x)

[j≠i , j≠i+1]

nulles sur cet

+

 La fonction Ni(x) est nulle partout sauf sur l’intervalle [xi-1 , xi+1] sur lequel elle est de forme

triangulaire. En généralisant, les fonctions Ni(x) sont de la forme :

x

x1 xi-1 xi xi+1  x1 x2 

x

x5 x6 

 N6(x) N1(x) Ni(x) , i=1 à 5

Fonctions "internes" Fonction à gauche Fonction à droite

Page 58: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 58/85

Modélisation

FN CRES

58

 

Ces fonctions sont appelées ‘Fonctions de Forme’. Elles permettent de reconstituer la

solution u(x) à partir de morceaux de droite sur chaque élément, connaissant les ui , valeur de

la fonction aux bornes de chaque élément.

 Remarque : en travaillant sur les ui , valeur de la fonction aux bornes de chaque élément, on

assure la continuité de la solution.

En choisissant les fonctions Ni(x) comme on vient de le voir, la solution s’exprime ainsi :

u x u N xi ii

( ) ( )==∑

1

6

 

La formulation variationnelle de l’équation à résoudre a établi que u(x) doit appartenir à

l’espace vectoriel V = { v : v∈L²(Ω) ; v’∈L²(Ω) ; v(0)=0; v(1)=0}.

La base des Ni(x) où on veut exprimer la solution est-elle un sous-espace vectoriel de V ?

Les Ni(x) et Ni

(x) appartiennent à L²(Ω), mais on remarque que tous les vecteurs ne sont pasnuls aux abscisses (0) et (1), en particulier N1 (x) et N6 (x).

On prendra donc comme sous-espace vectoriel pour l’approximation interne l’espace formé

 par la base de vecteurs N2 (x) à N5 (x), ce qui conduit a ne traiter que 4 inconnues. Ici, cette

conclusion est logique puisque les valeurs de la fonction en (0) et (1) sont imposées par des

conditions de Dirichlet homogènes et n’ont donc pas à faire partie des inconnues.

Ecriture algébrique

L’approximation de la solution dans l’espace vectoriel formé par les vecteurs Ni(x) retenus

conduit au système d’équations suivant (voir le paragraphe III-2.5).:

Trouver ui tel que ∀i =2,5 et ∀  j =2,5 : u a N N L Ni i ji

 j( , ) ( )=∑ =

2

5

 

avec a(Ni , N j) = ( ' ' )dx N N N Ni j i j+∫ 0

L(N j ) = f x N dx j( )0

1∫   

Calcul des termes a(Ni , N j)

Les fonctions Ni (x) étant identifiées, il n’y a pas de difficulté particulière à calculer ces

termes. Cependant, on peut remarquer certaines redondances dans les calculs.

• cas où | i - j | > 1 [i et j non contigus]

0

1

 Ni(x) N j(x)

x

xi-1 xi xi+1 x j-1 x j x j+1

 N(x)

 

On a [Ni (x) = 0 quand N j (x) ≠ 0] et [ N j (x) = 0 quand Ni (x) ≠ 0] d’où a(Ni , N j)=0.

Page 59: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 59/85

Modélisation

FN CRES

59

• cas où | i - j | = 1 [i et j contigus]

 Ni(x) N j(x)x

xi-1 xi xi+1=x j x j+1

 N(x)

1

0

 Il n’y a que l’intervalle [xi , x j ] où les fonctions sont non nulles simultanément et donc où

a(Ni , N j) est différent de zéro. Les fonctions étant identiques sur chaque intervalle, le calcul

mènera au même résultat sur chaque intervalle, en particulier sur un intervalle [0,h], h étant la

largeur d’un intervalle.

Soit sur h = xk+1 - xk , nous avons Nk (x) = 1−x

het Nk+1 (x) =

x

d’où a(Ni , N j) = ( ' ' ) [ ( ) ( )] N N N N dx

x

h

x

h h h dx

h

hi j i j

h h

+ = − + − = −∫ ∫ 0 0 1

1 1

6

1

 

• cas où | i - j | = 0 [i et j confondus]

1

 N j(x) Ni(x)

x

xi-1 xi=x j xi+1

 N(x)

 Il n’y a que l’intervalle [xi-1 , xi+1] où la fonction est non nulle et donc où a(Ni , Ni) est

différent de zéro. Cette intervalle correspond à une distance 2h sur laquelle la fonction Ni (x)

vautx

hsur [0,h] et 2−

x

hsur [h,2h].

d’où a(Ni , Ni) = ( ' ) ( ' ) [( ) ( ) ] N N dx N N dxx

h hdx

h

hi i

h

i i

h h2 2

0

2 2 2

0

2 2

02 2

1 2

3

2+ = + = + − = +∫ ∫ ∫   

AssemblageCette étape consiste à créer la matrice carrée du système matriciel en observant comment les

différents éléments sont agencés. On obtient :

u2 u3 u4 u5 

2

3

2h

h+   h

h6

1−   u2 

h

h6

1−   2

3

2h

h+   h

h6

1−   u3 

h

h6

1−   2

3

2h

h+   h

h6

1−   u4 

Page 60: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 60/85

Modélisation

FN CRES

60

 h

h6

1−   2

3

2h

h+   u5 

Calcul des termes L(Ni)

Connaissant les fonctions f(x) et Ni(x), il faut calculer L(Ni) = f x N x dxi( ) ( )0

1∫  .

Le calcul ne pose pas de difficulté particulière mais il peut être simplifier dans la mesure où

f(x) peut s'exprimer sous forme d'une fonction linéaire par morceaux.

Ici, en prenant une largeur h=0.2 sur notre exemple, nous avons :

2x - 1 = (-1) N1(x) + (-0.6) N2(x) + (-0.2) N3(x) + (0.2) N4(x) + (0.6) N5(x) + (1) N6(x)

soit f(x) = βk k 

 N x( )

=

∑1

6

 

d'où L(Ni) = [ ( )] ( )βk k k 

i N x N x dx=

∑∫ 1

6

0

1= [ ( ) ( ) ]βk k i

 N x N x dx0

1

1

6

∫ ∑=

 

or les expressions N x N x dxk i( ) ( ) ont déjà été calculées auparavant sur un élément et par 

ailleurs, Ni (x) est nulle en dehors de l’intervalle [xi-1 , xi+1]

On a alors :

L(Ni) = βi-1 N x N x dxi i

h

−∫  10( ) ( ) + βi N x N x dxi i

h( ) ( )

0

2∫  + βi+1 N x N x dxi i

h

+∫  10( ) ( )

L(Ni) = βi-1 ( h6

) + βi ( 23h ) + βi+1 ( h

6)

L’application numérique donne :

L(N )= -0.12 ; L(N )= -0.04 ; L(N )= 0.04 ; L(N )= 0.12 ;

Résolution du système

Elle se fait pour h=0.2

6.5.2. Cas de condition de Neumann homogèneOn reprend l'exemple précédent et les conditions limites suivantes :

-u"(x) + u(x) = f(x) = 2x - 1 sur le domaine Ω = [0,1]

avec u(0) = 0

u’(1) = 0

On a vu que ce problème a la même solution que le problème suivant :

trouver u∈V tel que ∀ v∈V , on ait : a(u,v) = L(v)

Page 61: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 61/85

Modélisation

FN CRES

61

avec V = { v : v∈L²(Ω) ; v’∈L²(Ω) ; v(0)=0}

a(u,v) = ( ' ' )uv u v dx+∫ 01

 

L(v) = ( )fv dx0

1

∫   

Les développements sont les mêmes qu’au paragraphe précédent. Ici, la base pour 

l’approximation interne sera constituée des vecteurs N2(x) à N6(x) car N6(x)∈V [la condition

v(1)=0 d’appartenance à V ayant disparue]

On prendra simplement garde au fait que a(N6, N6) ne s’intègre que sur un seul élément car la

fonction N6 est différente de zéro que sur le dernier intervalle. On aura donc a(N6, N6)=h

h3

1+  

et la matrice carrée s’écrit :

u2 u3 u4 u5 u6 

23

2hh

+   hh61−   u2 

h

h6

1−   2

3

2h

h+   h

h6

1−   u3 

h

h6

1−   2

3

2h

h+   h

h6

1−   u4 

h

h6

1

−  2

3

2h

h+  h

h6

1

−   u5

h

h6

1−   h

h3

1+   u6 

Le calcul des L(Ni) est identique, on rajoute le calcul de L(N6) qui donne 0.0866.

6.5.3. Cas de condition de Neumann non homogèneOn reprend l'exemple précédent et les conditions limites suivantes :

-u"(x) + u(x) = f(x) = 2x - 1 sur le domaine Ω = [0,1]

avec u(0) = 0

u’(1) = g N 

On a vu que ce problème a la même solution que le problème suivant :

trouver u∈V tel que ∀ v∈V , on ait : a(u,v) = L(v)

avec V = { v : v∈L²(Ω) ; v’∈L²(Ω) ; v(0)=0}

a(u,v) = ( ' ' )uv u v dx+∫ 01

 

Page 62: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 62/85

Modélisation

FN CRES

62

L(v) = ( )fv dx0

1

∫  + v(1)g N 

Les termes a(Ni ,N j ) sont identiques aux précédents.

Les termes L(Ni ) s’écrivent : L(Ni) = f x N x dxi( ) ( )0

1∫  + Ni (1)g N.

Le dernier terme est toujours nul sauf pour i=6 où N6(1)=1. Aussi, les L(Ni ) ne changent pas

non plus sauf L(N6 ) auquel on rajoute la quantité g N.

6.5.4. Cas de condition de Dirichlet non homogèneOn reprend l'exemple précédent et les conditions limites suivantes :

-u"(x) + u(x) = f(x) = 2x - 1 sur le domaine Ω = [0,1]

avec u(0) = 0

u(1) = gD 

On a vu que ce problème a la même solution que le problème suivant :

Soit u = u0 + w la solution de l’équation ; u0∈U0 

trouver w∈V tel que ∀v∈V, on ait a(w,v) = L(v)

avec U0 = { u0 : u0∈L²(a,b); u0’∈L²(a,b); u0(a)=gD }

V = { v : v∈L²(a,b); v’∈L²(a,b); v(a)=0 }

a(w,v) = ( 'va

 bw' + vw)dx∫  et L(v) = ( )dxvf 

a

 b

∫   

En fait, on peut considérer que les deux espaces U0 et V sont engendrés par la même base de

vecteurs Ni (x) et sont donc de la forme αi ii

n

 N x( )=∑1

, mais pour 

- U0 : α1 = gD et αi = 0 pour les autres coefficients

- V : α1 = u1 =0 et αi = ui ≠ 0 (à priori) pour les autres coefficients

Chacun des vecteurs u0 [∈ U0] et w [∈ V] peut alors s’écrire sur cette base.

Un vecteur solution u = u0 [∈ U0] + w [∈ V] s’écrit alors :

u = βi ii

 N x( )=∑

1

5

[∈ U0] + u N xi ii

( )=∑

1

5

[∈ V]

Le système d’équations (voir le paragraphe III-2.5) auquel conduit l’approximation interne va

alors s’écrire :

∀  j =1,n : βi i ji

i i ji

 ja N N u a N N L N( , ) ( , ) ( )= =∑ ∑+ =

1

5

1

5

 

Or pour U0 : β1 = gD et βi = 0 pour les autres coefficients

  pour V : u1 = 0 et ui ≠ 0 (à priori) pour les autres coefficients

La recherche des inconnus ui conduit à un problème à 4 inconnues pour i=2,3,4 et 5 soit :

∀ j =2,n : g a N N u a N N L ND i i ji

 j( , ) ( , ) ( )1 22

5

+ ==∑  

∀  j =2,n : u a N N L N g a N Ni i ji

 j D( , ) ( ) ( , )=∑ = −2

5

1 2  

Page 63: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 63/85

Modélisation

FN CRES

63

Par rapport au système d’équations du premier exemple (paragraphe III-4.1.1) seule la

 première ligne est changée et s’écrit :

u2 a(N2 ,N2 ) + u3 a(N2 ,N3 ) + 0 = L(N2 ) - gD a(N1 ,N2 )

6.6. Cas d’une équation différentielle - Fonction parabolique parmorceaux

On reprend l'exemple précédent avec les conditions limites suivantes :

-u"(x) + u(x) = f(x) = 2x - 1

avec u(0) = 0

u'(1) = 0

On a vu que ce problème a la même solution que le problème suivant :

trouver u∈V tel que ∀ v∈V , on ait : a(u,v) = L(v)

avec V = { v : v∈L²(Ω) ; v’∈L²(Ω) ; v(0)=0}

a(u,v) = ( ' ' )uv u v dx+∫ 01

 

L(v) = ( )fv dx0

1

∫   

Choix des éléments

Le domaine d’étude est séparé en éléments de largeur (h). Cette largeur peut être différente

 pour chaque élément (dans le cas où on veut affiner la connaissance de la solution à certains

endroits)

Dans cet exemple, on choisit de séparer le domaine d'étude en 5 éléments de même longueur.

x1 x2 x3 x4 x5 x6

 

Approximation parabolique par morceaux

On approxime la solution sur chaque élément par un morceau de parabole

Page 64: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 64/85

Modélisation

FN CRES

64

 

x 1 x 2 x 3 x 4 x 5 x 6 

 Remarque :

- il y a continuité de la fonction par morceaux

- il n'y a pas continuité de la dérivée (comme pour l'approximation linéaire par 

morceaux)

Pour décrire une parabole sur un élément, il suffit de 3 points; on choisit généralement les

deux points extrêmes (ce qui assurera la continuité entre les morceaux) et le point milieu.

xi xi+1/2 xi+1

 La dimension de l'espace vectoriel où se fait l'approximation est déterminée par le nombre

d'inconnues que l'on se fixe. Ici, 3 inconnues par éléments, dont les inconnues extrêmes sont

communes à deux autres éléments (aux conditions limites près que l'on verra un peu plusloin).

La solution u(x) sera donc approximée par l'expression suivante exprimée sur une base de 11

vecteurs -à priori- notés Ni(x) et Ni+1/2(x)

u(x) = u N xi

i

i

=∑

1

5

( ) + u N xi

i

i+=

+∑ 1 2

1

5

1 2/ / ( ) 

Le premier terme de la somme correspond aux points extrêmes de chaque élément; le

deuxième terme correspond aux points milieux.

Le problème est alors le choix des fonctions Ni(x) et Ni+1/2(x)

En adoptant le même symbolisme vue avec le paragraphe précédent, on peut écrire :

Page 65: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 65/85

Modélisation

FN CRES

65

 

= ui  + ui+1/2 + ui+1 

xi xi+1/2 xi+1 xi xi+1/2 xi+1 xi xi+1/2 xi+1 xi xi+1/2 xi+1

u(x)  Ni(x)  Ni+1/2(x) Ni+1(x)

1 1 1

 

Les fonctions Ni sont donc de la forme :

x1 x1/2 x2 xi-1 xi-1/2 xi xi+1/2 xi+1 xi xi+1/2 xi+1

 N1 Ni   N6 

1

Les fonctions Ni+1/2 sont de la forme :

xi xi+1/2 xi+1  

 Ni+1/2 

1

Sur 1 élément donné, il y a 3 fonctions qui sont simultanément non nulles :

xi xi+1/2 xi+1 

 Ni+1/2

 Ni   Ni+1 

Page 66: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 66/85

Modélisation

FN CRES

66

On remarque que sur la frontière x1 où on a une condition de Dirichlet, le vecteur N1 

n'appartient pas à l'espace vectoriel où on approxime la solution.

En effet, cet espace est tel que ∀ v : v(0) = 0; or N1(0) = 1;

L'espace vectoriel sera donc constitué à partir de la base :

( N2 , N3 , N4 , N5 , N6 , N3/2 , N5/2 , N7/2 , N9/2, N11/2 )

Ecriture algébrique

En remplaçant u(x) et v(x) par leurs expressions sur la base d'approximation :

u(x) = u N xi

i

i

=∑

1

5

( ) + u N xi

i

i+=

+∑ 1 2

1

5

1 2/ / ( ) 

v(x) = α i

i

i N x

=

∑1

5

( ) + α i

i

i N x+

=

+∑ 1 2

1

5

1 2/ / ( ) 

et en utilisant la linéarité de la forme L , l'équation a(u,v) = L(v) se ramène au système

d'équations suivant avec 10 inconnues :

(en notant Ni le vecteur Ni(x) et Ni+1/2 le vecteur Ni+1/2(x) )

∀ i = 2, 6 a ( u N j

 j

 j

=∑

2

6

+ u N j

 j

 j+=

+∑ 1 2

1

5

1 2/ / , Ni ) = L(Ni)

∀ i = 1, 5 a ( u N j

 j

 j

=∑

2

6

+ u N j

 j

 j+=

+∑ 1 2

1

5

1 2/ / , Ni+1/2 ) = L(Ni+1/2)

soit encore en utilisant la bilinéarité de la forme notée a :

∀ i = 2, 6 u a N N j

 j

 j i

=∑

2

6

( , ) + u a N N j

 j

 j i+=

+∑ 1 2

1

5

1 2/ /( , ) = L(Ni)

∀ i = 1, 5 u a N N j

 j

 j i

=+∑

2

6

1 2( , )/ + u a N N j

 j

 j i+=

+ +∑ 1 2

1

5

1 2 1 2/ / /( , ) = L(Ni+1/2)

Ce système d’équations correspond au système matriciel suivant :

a(N N a(N N a(N N a(N N

a(N N a(N N a(N N a(N N

a(N N a(N N a(N N a(N N

a(N N a(N N a(N N

3 2 3 2 3 2 2 3 2 5 2 3 2 6

2 3 2 2 2 2 5 2 2 6

5 2 3 2 5 2 2 5 2 5 2 5 2 6

11 2 3 2 11 2 2 11 2 5

/ / / / / /

/ /

/ / / / / /

/ / / / /

, ) , ) , ) ... , )

, ) , ) , ) ... , )

, ) , ) , ) ... , )

... ...

... ...

... ...

... ...

... ...

, ) , ) , 2 11 2 6

6 3 2 6 2 6 5 2 6 6

) ... , )

, ) , ) , ) ... , )

/

/ /

a(N N

a(N N a(N N a(N N a(N N

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

u

u

u

u

u

u

u

u

u

u

3 2

2

5 2

3

7 2

4

9 2

5

11 2

6

/

/

/

/

/

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

=

L N

L N

L N

L N

L N

( )

( )

( )

...

...

...

...

...

( )

( )

/

/

/

3 2

2

5 2

11 2

6

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

 

Page 67: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 67/85

Modélisation

FN CRES

67

 

Calcul des termes a(Nk,Nm)

Sachant que a est symétrique, l'écriture de ce système impose de calculer les expressions de

la forme suivante :

a(Ni,N j) = ( ' ' ) N N N N dxi j i j+∫ 01

 

a(Ni+1/2,N j) = ( ' ' )/ / N N N N dxi j i j+ ++∫  1 2 1 20

1

 

a(Ni+1/2,N j+1/2) = ( ' ' )/ / / / N N N N dxi j i j+ + + ++∫  1 2 1 2 1 2 1 20

1

 

Chaque calcul d'intégrale des expressions ci-dessus entre 0 et 1 peut être décomposé sur 

chaque élément; en effet (en notant g(x) un des termes à intégrer ci-dessus) :

g x dx( )0

1

∫  = g xx

x

( )dx1

2

∫  + g xx

x

( )dx2

3

∫  + ... + g xx

x

( )dx5

6

∫  = g xx

x

i i

i

( )dx+

∫ ∑=

1

1

5

 

soit g x dx( )0

1

∫  = g xh

( )dx0∫  + g x

h

h

( )dx2

∫  + ... + g xh

h

( )dx4

5

∫  = g xe

é lé ment

( )dx( )∫ ∑  

Comme toutes les fonctions Ni et Ni+1/2 ont la même allure quelque soit l'élément où on se

  place, les calculs d'intégrales sont les mêmes sur chaque élément et en particulier sur le

 premier entre 0 et h (de plus, on a pris dans cet exemple tous les éléments de même largeur h).

On calcule donc chaque terme intégral entre 0 et h et on procédera ensuite à l'assemblage

des éléments.

Calcul des équations des morceaux de parabole

L'équation générale s'écrit : ax2 + bx + c = 0. On sait que chaque parabole passe par 3 points

ce qui donne 3 équations à 3 inconnues à chaque fois.

0  h/2  h

 Ni  passe par (0,1), (

h

2, 0) et (h,0)

soit l'équation Ni(x) =2

h² x² -3

h x + 1 = 0

0 h/2 h

 Ni+1/2

 

 passe par (0,0), (h

2,1) et (h,0)

soit l'équation Ni+1/2(x) = −4

h²x² +

4

hx = 0

Page 68: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 68/85

Modélisation

FN CRES

68

 

0 h/2 h

 Ni+1   passe par (0,0), (h

2,0) et (h,1)

soit l'équation Ni+1(x) =2

h²x² -

1

hx = 0

Maintenant que l'on connaît les équations des branches de paraboles, il faut calculer les

termes a(Nk ,Nm), ce qui impose de calculer des intégrales avec les fonctions du type Nk et N'k .

Les calculs sont récapitulés dans les tableaux suivants :

 Ni Ni+1/2 Ni+1 

( ) N N dxk m

h

0∫    Ni  2

15

h  h

15  −

h

30 

 Ni+1/2  h

15  8

15

h   h

15 

 Ni+1 −h

30  h

15  2

15

 N'i N'i+1/2 N'i+1 

( ' ' ) N N dxk m

h

0∫    N'i  7

3h

  −8

3h

  1

3h

 

 N'i+1/2  −8

3h  16

3h  −

8

3h 

 N'i+1  1

3h  −

8

3h  7

3h 

Assemblage

 fonction N i+1/2 

Chaque fonction Ni+1/2 est non nulle sur un seul élément à la fois. Sur cet élément, deux autresfonctions, Ni et Ni+1 sont également non nulles

L'inconnue rattachée à la fonction en question, ui+1/2 , sera donc liée par une équation avec lesinconnues ui et ui+1; Il n'y aura donc que 3 termes sur une ligne de la matrice relative à

l'équation (i+1/2) d'un point milieu.

 fonction N i 

Chaque fonction Ni est non nulle sur deux éléments à la fois. Sur ces éléments, quatre autresfonctions, Ni-1, Ni-1/2, Ni+1/2 et Ni+1 sont également non nulles

Page 69: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 69/85

Modélisation

FN CRES

69

L'inconnue rattachée à la fonction en question, ui , sera donc liée par une équation avec les

inconnues ui-1, ui-1/2, ui+1/2, et ui+1; Il y aura donc 5 termes sur une ligne de la matrice relative àl'équation (i) d'un point extrême d'élément.

Les points limites à droite et à gauche constituent des cas particuliers.En x3/2 , il n'y aura que deux termes dans l'équation car N1(x) n'existe pas dans la base

d'approximation.

En x6 , qui est un point extrême, mais frontière, il n'y a pas d'élément à droite. Il y aura donc 3

termes à considérer (relatifs à N5, N11/2 et N6).

 Remarque :

Attention aux points extrêmes non frontières pour lesquels deux éléments sont concernés; ilfaut cumuler deux fois les intégrales du type a(Ni,Ni) [puisque chaque calcul n'a été fait que

sur un seul élément]

On en déduit les coefficients de la matrice carrée du système à résoudre :

u3/2 u2 u5/2 u3 u7/2 u4 u9/2 u5 u11/2 u6 

8

15

16

3

h

h+   h

h15

8

3−   u3/2 

h

h15

8

3−   4

15

14

3

h

h+   h

h15

8

3−   − +

h

h30

1

3

 

u2 

h

h15

8

3−  8

15

16

3

h

h+  h

h15

8

3− 

u5/2 

− +h

h30

1

3

 

h

h15

8

3−   4

15

14

3

h

h+   h

h15

8

3−

 

− +h

h30

1

3

 

u3 

h

h15

8

3−   8

15

16

3

h

h+

 

h

h15

8

3−

 

u7/2 

− +h

h30

1

h

h15

8

3−

 

4

15

14

3

h

h+

 

h

h15

8

3−

 

− +h

h30

1

u4 

h

h15

8

3−

 

8

15

16

3

h

h+

 

h

h15

8

3−

 

u9/2 

− +h

h30

1

3

 

h

h15

8

3−

 

4

15

14

3

h

h+

 

h

h15

8

3−   − +

h

h30

1

3

 

u5 

h

h15

8

3

 

8

15

16

3

h

h

+   h

h15

8

3

−   u11/2 

Page 70: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 70/85

Modélisation

FN CRES

70

− +h

h30

1

3

 

h

h15

8

3−   2

15

7

3

h

h+   u6 

Calcul des termes L(Ni) et L(N

i+1/2)

Connaissant les fonctions f(x) et Nk (x), il faut calculer L(Nk ) = f x N x dxk ( ) ( )0

1

∫  .

Ici, le calcul peut être simplifier dans la mesure où f(x) peut s'exprimer sous forme d'unefonction parabolique par morceaux.

En effet :2x - 1 = (-1)N1(x) + (-0.8)N3/2(x) + (-0.6)N2(x) + (-0.4)N5/2(x) + (-0.2)N3(x) + (0)N7/2(x) +

(0.2)N4(x) + (0.4)N9/2(x) + (0.6)N5(x) + (0.8)N11/2(x) + (1)N6(x)

soit f(x) = βm m

m

N x( )=

∑1

11

 

d'où L(Nk ) = [ ( )] ( )βm m

m

k  N x N x dx=∑∫ 1

11

0

1 = [ ( ) ( ) ]βm m k 

m

 N x N x dx0

1

1

11

∫ ∑=

 

or les expressions N x N x dxm kxi

xi

( ) ( )+∫  1

ont déjà été calculées auparavant sur un élément.

On a par exemple :

L(N2) = (-1) N x N x dxh

1 20

( ) ( )∫  + (-0.8) N x N x dxh

3 2 20

/ ( ) ( )∫  + (-0.6) N x N x dxh

2 20

2

( ) ( )∫  +

(-0.4) N x N x dxh

5 2 20

/ ( ) ( )∫  + (-0.2) N x N x dxh

3 20

( ) ( )∫   

(On ne fait pas intervenir les autres fonctions car N2(x) est alors égal à zéro.)

L(N2) = (-1)(−h

30) + (-0.8)(

h

15) + 2(-0.6)(

2

15

h) + (-0.4)(

h

15) + (-0.2)(−

h

30) = −

h

On obtient

L(N3/2) = −8

15

hL(N2) = −

h

5L(N5/2) = −

4

15

hL(N3) = −

h

15 

L(N7/2) = 0 L(N4) =h

15L(N9/2) =

4

15

hL(N5) =

h

L(N11/2) =8

15

hL(N6) =

h

Résolution du système

Elle se fait pour h=0.2

6.7. Cas d’une EDP – Approximation linéaire par élément

Soit le problème :

Page 71: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 71/85

Modélisation

FN CRES

71

 

u=0

u=0 u=0

u=0

-Δu = f 

100

1

x

y

soit Ω = (]0,1[ x ]0,1[) le domaine interne et Γ la frontière.On montre en utilisant la formule de Green :

( . ) ( ) ( )grad u grad v d v u d vu

nd

 ⎯ →  ⎯ ⎯ → ⎯ 

∫ ∫ ∫ + =Ω Ω Γ

Ω Δ Ω Γ∂∂

et en choisisant v/Γ = 0

que l'équation (− =Δu f ) admet la même solution que le problème suivant :

( . ) ( )grad u grad v d v f d ⎯ →  ⎯ ⎯ → ⎯ 

∫ ∫ =Ω ΩΩ Ω

 

ce qui peut aussi s'écrire :

( ) ( )∂∂

∂∂

∂∂

∂∂

u

x

v

x

u

y

v

ydx dy v f dx dyΩ Ω∫∫ ∫∫  + =  

d'où le nouveau problème : trouver u∈V tel que ∀v∈V , on ait a(u,v) = L(v)

avec V = {v : v∈L²(Ω) ;∂∂v

x∈L²(Ω) ;

∂∂v

y∈L²(Ω) ; v/Γ=0} = H 0

1 (Ω)

a(u,v) = ( . )grad u grad v d ⎯ →  ⎯ ⎯ → ⎯ 

∫  ΩΩ

 

L(v) = ( )v f dΩΩ∫   

Choix des éléments

Parmi les plus simples à mettre en oeuvre, il y a les éléments triangulaires (c’est une surface

descriptible à partir d’un nombre de points minimum : 3).

Page 72: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 72/85

Modélisation

FN CRES

72

100

1

1 2 3 4

9 10 11 12

5 6 7 8

13 14 15 16

 

Il y a en tout dans cet exemple 50 éléments triangulaires qui permettre de couvrir tout ledomaine d’étude.

(Sur ce schéma, les noeuds qui vont nous intéresser sont numéroter de 1 à 16.)

Approximation de la solution

On veut exprimer la solution dans un sous-espace vectoriel de vecteurs de base Ni(x) tels que

u(x) = u N x yi i

i

m

( , )=

∑1

avec Ni(x,y) fonctions linéaires par morceaux en x et en y

Il suffit de choisir les Ni(x,y) de la forme suivante :

1

i

 Ni(x)x

y

  Ni(x,y) vaut (1) au sommet (i) et (0) sur les éléments triangulaires qui ne sont pas jointifs avec

le sommet (i).

On obtient aux bords :

1

iiix

y

 

Page 73: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 73/85

Modélisation

FN CRES

73

Les vecteurs sur la frontière du problème n'appartiennent pas à ce sous-espace vectoriel car 

 Ni(x,y) = 1 sur Γ alors qu'on doit avoir = 0 sur Γ 

On choisit donc comme base l'ensemble des vecteurs définis sur les noeuds réguliers

(internes), soit ici 16 vecteurs correspondants aux noeuds numérotés de 1 à 16 sur le schéma

 précédent.

Ecriture algébrique

En utilisant la linéarité de a(u,v) et de L(v), on se ramène au système d'équations suivant :

trouver u j avec j = 1.. 16 tq ∀i=1..16 : u a N N L N j i j j

i( , ) ( )==

∑1

16

 

Calcul des termes a(Ni,N j)

L'intégrale sur Ω peut se décomposer sur chaque élément :

a(Ni,N j) = grad N grad N dxdyi j

 ⎯ →  ⎯ ⎯ → ⎯ 

∫∫  .Ω

= ( . )( )

grad N grad N dxdyi je

é lé ment

 ⎯ →  ⎯ ⎯ → ⎯ 

∫∫ ∑  

On va donc déterminer les intégrales sur chaque triangle (ce sont les mêmes) et ensuite faire

l'assemblage.

Sur chaque triangle, il y a 3 fonctions qui sont simultanément non nulles. De plus, il y a 2

types de triangle :

type G

type D

 b

a

ch

-h

 b

c

a

h

-h

 

Triangle de type Gil faut d'abord calculer l'équation des fonctions non nulles sur ce triangle. Calculons par 

exemple la fonction notée N b(x,y) valant 1 au sommet (b).

 N b(x,y) = b1 + b2 x + b3 y avec N b(0,0) = 1

 N b(h,0) = 0

 N b(0,-h) = 0

On en tire N b(x,y) = 1 -x

h+

y

h+ d'où grad

 ⎯ → ⎯ 

N b =−⎡

⎣⎢

⎦⎥

1

1

h

De façon analogue

 Na(x,y) = - yh

d'où grad ⎯ → ⎯  Na = 01−

⎡⎣⎢ ⎤

⎦⎥h 

Page 74: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 74/85

Modélisation

FN CRES

74

 

 Nc(x,y) =x

hd'où grad

 ⎯ → ⎯ 

Nc =1

0

h⎡

⎣⎢

⎦⎥ 

d'où le calcul des intégrales dans ce triangle (on a y = x-h)

exemple : a(Na, Na) = ( )2 2

0

0

 Na

x

 Na

ydy dx

x h

h

+⎡

⎢⎢

⎥⎥

⎛ ⎝ ⎜ ⎞

 ⎠⎟∫ ∫  −

= (²)0

10

0+⎡

⎣⎢⎤⎦⎥−∫ ∫  h

dy dxx h

h

 

= (²

)y

hdx

h

x h

⎡⎣⎢

⎤⎦⎥∫ 

−0

0

=h x

hdx

h −∫  ²0

=x

h

x

h

h

−⎡⎣⎢

⎤⎦⎥

²

²2 0

= 1 -1

2=

1

Tous calculs faits :

sommets a b c

a1

2  −

1

2   0

 b −1

2  1 −

1

c  0 −1

2  1

triangle de type G

Triangle de type D

De façon similaire

grad ⎯ → ⎯ 

Na =−⎡

⎣⎢

⎦⎥

1

0

hgrad ⎯ → ⎯ 

N b =1

1

h

h−⎡

⎣⎢

⎦⎥ grad

 ⎯ → ⎯ 

Nc =0

1 h

⎣⎢

⎦⎥ 

sommets a b c

a 1

2

  −1

2

  0

 b −1

2  1 −

1

c  0 −1

2  1

triangle de type D

Assemblage

Page 75: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 75/85

Modélisation

FN CRES

75

 

Il s'agit de représenter le terme ( . )( )

grad N grad N dxdyi je

é lément

 ⎯ →  ⎯ ⎯ → ⎯ 

∫∫ ∑ et donc de cumuler les

intégrales sur chaque élément (il y a ici 50 triangles).

Cette opération se fait en observant comment s'organisent les triangles.

Calcul de a(N i , N  j) pour i=jOn adopte la numérotation suivante des éléments quand on se trouve sur le sommet (i) :

i

1

2

3

4

5

6

 

Seuls ces 6 éléments interviennent puisque la fonction N i(x,y) est non nulle sur ceux-ci, et est

nulle sur les autres. Il faut donc regarder le type (G ou D) de triangle concerné à chaque fois.

élément ou triangle

type sur ce triangle, (i) est

un sommet de type

valeur de a(Ni,Ni) sur 

cet élément

1 D c 1/2

2 G b 1

3 D a 1/24 G a 1/2

5 D b 1

6 G c 1/2

d’où a(Ni,Ni) = 4

Calcul de a(N i , Nj) pour | i-j | = 1

oui

1

2

i+1

cas A

i

i+1

2

1

cas C

oui+1

1

2

i

cas A’Frontièregauche

Frontière

droite  

Remarque : Il y a aussi un cas C’ a envisager avec (i+1) en bas et (i) en haut. Mais

a(Ni, Nj) étant symétrique, on trouvera la même chose dans les cas A et A’, et C et

C’ On ne calcule que le cas A et C

Page 76: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 76/85

Modélisation

FN CRES

76

Cas A .

élément ou triangle

type (i) et (i+1) sont des

sommets de type

valeur de a(Ni,Ni) sur 

cet élément

cas A 1 G b et c -1/2

cas A 2 D a et b -1/2

d’où a(Ni,Ni+1) = a(Ni+1,Ni) = -1 dans les cas A et A’

Cas C

Les fonctions Ni et Ni+1 ne sont jamais non nulles simultanément,

d’où a(Ni,Ni+1) = a(Ni+1,Ni) = 0 dans ce cas

Calcul de a(N i , Nj) pour | i-j | = 4 (en généralisant, =n si il y a n sommets par ligne)

i

1

2

i+1

ou

i+1

1

2

icas A cas A’

 Encore une fois, les deux cas sont symétriques

élément ou triangle

type (i) et (i+1) sont des

sommets de type

valeur de a(Ni,Ni) sur 

cet élément

cas A 1 D b et c -1/2

cas A 2 G a et b -1/2

d’où a(Ni,Ni+1) = a(Ni+1,Ni) = -1 dans ces cas.

Calcul de a(N i , Nj) pour | i-j | = 4+1 = 5 (en généralisant,=n+1 si il y a n sommets par ligne)

ou

cas A cas A’

i

1

2

i+1

i+1

1

2

i

 Les deux cas sont toujours symétriques

élément ou triangle

type (i) et (i+1) sont des

sommets de type

valeur de a(Ni,Ni) sur 

cet élément

cas A 1 D a et c 0

cas A 2 G a et c 0

d’où a(Ni,Ni+1) = a(Ni+1,Ni) = 0 dans ces cas.

Calcul de a(N i , N  j ) pour tous les autres cas

Page 77: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 77/85

Modélisation

FN CRES

77

Il y a toujours au moins une des fonctions qui est nulle d’où a(Ni,N j) = 0 pour ces cas là.

On peut récapituler ces résultats dans une matrice.

u1 u2 u3 u4 u5 u6 u7 u8 u9 u10 u11 u12 u13 u14 u15 u16 4 -1 -1 u1 

-1 4 -1 -1 u2 

-1 4 -1 -1 u3 

-1 4 0 -1 u4 

-1 0 4 -1 -1 u5 

-1 -1 4 -1 -1 u6 

-1 -1 4 -1 -1 u7 

-1 -1 4 0 -1 u8 

-1 0 4 -1 -1 u9 

-1 -1 4 -1 -1 u10 -1 -1 4 -1 -1 u11 

-1 -1 4 0 -1 u12 

-1 0 4 -1 u13 

-1 -1 4 -1 u14 

-1 -1 4 -1 u15 

-1 -1 4 u16 

Calcul des termes L(Ni)

Ces termes s’écrivent L(Ni) = f x N x dx dyi( ) ( )0

1

0

1

∫ ∫  . Leurs calculs ne posent pas de difficulté

 particulière, connaissant f(x) et les fonctions Ni(x).

6.8. Cas d’une EDP – Approximation quadratique par élément

 Nous indiquons juste ici les fonctions de forme dans le cas d’éléments triangulaires. Il y en a

deux. Les points noirs indiquent les valeurs remarquables.

Page 78: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 78/85

Modélisation

FN CRES

78

1

1

6.9. Extension de la méthode

A partir des principes évoqués, on peut mettre en œuvre des fonctions de formes de type

cubiques, voire quartiques. Les calculs deviennent alors très fastidieux et sortent du cadre de

ce cours qui est une initiation à cette méthode.

De même, on peut traiter ces problèmes avec des éléments rectangulaires.

Pour les problèmes à trois dimensions d’espaces, la méthode s’étend sans difficulté : les

triangles sont remplacés par des tétraèdres et les rectangles par des parallélépipèdes.

6.10. Extension aux problèmes évolutifs

L’introduction du temps dans l’équation se traite comme pour la méthode des différences

finies, avec des schémas de type explicites et implicites.

Page 79: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 79/85

Modélisation

FN CRES

79

 

7. RESOLUTION DES SYSTEMES D’EQUATIONS LINEAIRESCe chapitre est là à titre informatif, car il permet de compléter ce cours sur la modélisation.

7.1. IntroductionOn cherche à résoudre le système suivant :

⎪⎪⎩

⎪⎪⎨

=+++

=+++

=+++

nnnn22n11n

2nn2222121

1nn1212111

 bxaxaxa

 bxaxaxa

 bxaxaxa

L

M

L

L

 

On note ce système

BAX = où ijaA = , [ ]ixX = et [ ]i bB =  

Soit

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

n

2

1

n

2

1

nn2n1n

n22221

n11211

 b

 b

 b

x

x

x

 

aaa

aaa

aaa

MM

L

M

L

L

 

On ne change rien :

-  en permutant deux lignes simultanément en A et B

-  en multipliant chaque ligne par un scalaire

-  en additionnant à une ligne une autre ligne de la matrice-  en combinant les opérations précédentes

 Remarque : il existe une seule solution X si la matrice A est régulière (déterminant non nul)

La nécessité de résoudre de tels systèmes d'équations apparaît dans de nombreuses méthodes

numériques (analyses de données, résolutions d'équations, méthodes des différences finies et

des éléments finis,…).

On peut distinguer deux principes de résolution d'un tel système :

-  la méthode directe (basée sur la méthode du pivot) qui s'adapte bien aux systèmes

où A est une matrice dense mais d'ordre peu élevé, ou bien A est une matrice bande, même d'ordre élevé.

-  La méthode itérative, qui s'adapte bien aux systèmes où A est creuse (beaucoup de

zéros) et d'ordre élevé.

7.2. Méthodes directes  Nous allons traiter ici la méthode du pivot et ses variantes. Le but est de transformer le

système en un système résolvable directement.

Page 80: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 80/85

Modélisation

FN CRES

80

7.2.1. Méthode de Gauss-Jordan ou méthode du pivot7.2.1.1  Principe

Soit le système

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

n

2

1

n

2

1

nn2n1n

n22221

n11211

 b

 b

 b

x

x

x

 

aaa

aaa

aaa

MM

L

M

L

L

 

on cherche a le transformer en un système de la forme

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

n

2

1

n

2

1

d

d

d

x

x

x

 

100

010

001

MM

L

M

L

L

 

On cherche donc à transformer la matrice A en la matrice identité.

a)  on divise la première ligne par  11a  

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

n

2

1

n

2

1

nn2n1n

n22221

n112

 b

 b

' b

x

x

x

 

aaa

aaa

'a'a1

MM

L

M

L

L

avec

⎪⎪⎩

⎪⎪⎨

=

=

11

11

11

 j1

 j1

a

 b' b

a

a'a

pour j = 1 à n

 b)  on annule le premier terme des autres lignes en retranchant à chaque ligne la première

multipliée par le premier élément de chaque ligne.

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

n

2

1

n

2

1

nn2n

n222

n112

'' b

'' b

' b

x

x

x

 

''a''a0

''a''a0

'a'a1

MM

L

M

L

L

avec⎩⎨⎧

−=

−=

1i1ii

1i j1ijij

a ' b b'' b

a 'aa''apour i = 2 à n et j = 1 à n

On dit que l'on vient de prendre la première ligne comme pivot et son premier élément

comme élément pivot.On recommence en prenant comme pivot la seconde ligne et son deuxième élément comme

élément pivot.

On arrive ainsi au système

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

n

2

1

n

2

1

d

d

d

x

x

x

 

100

010

001

MM

L

M

L

L

 

d'où la résolution de l'équation avec

ii dx = pour i = 1 à n

Page 81: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 81/85

Modélisation

FN CRES

81

7.2.1.2   Limite de la méthodeSi l'élément pivot est nul (ou proche de 0), on ne peut pas faire la division (le résultat peut être

faussé).

Solutions pour la programmation:

a)  Solution simpleOn teste la valeur du pivot et on n'effectue l'opération que si elle est supérieure à une valeur  ε,

sinon on arrête le calcul.

b)   Recherche du pivot maximum partielParmi les lignes pivot qui n'ont pas été traitées, on cherche la ligne (i) pour laquelle l'élément

 pivot est maximum. Puis on intervertit les lignes pour travailler avec le pivot maximum

(en valeur absolue), sans oublier d'intervertir les termes ib .

Ou bien on procède de la même façon en intervertissant les colonnes de façon à placer en

élément pivot le terme maximum de la ligne pivot (sans oublier d'inverser les ix ).

c)   Recherche du pivot totalOn cherche la valeur maximale de la matrice (*) pour s'en servir d'élément pivot, en

intervertissant les lignes et les colonnes.

(*) il s’agit ici de la sous-matrice carrée constituée de l’intersection des lignes et colonnes

non encore traitées. 

Remarque : si le dernier pivot est nul, alors c'est que le système d'équation admet une infinité

de solutions.

d)  Vérification des résultats

On peut aussi vérifier les résultats en recalculant les équations (connaissant les ix ) et les

comparer aux termes ib .

7.2.2. Méthode de Gauss (ou triangularisation)On cherche cette fois à se ramener à un système de la forme :

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

n

2

1

n

2

1

n2

n112

d

d

d

x

x

x

 

100

c10

cc1

MM

L

M

L

L

 

On cherche donc à transformer la matrice A en une matrice triangulaire supérieure

équivalente (on peut aussi raisonner avec une matrice triangulaire inférieure).

L'algorithme est le même que pour la méthode précédente, mais l'étape de soustraction ne

s'effectue que pour les lignes inférieures.

Puis on trouve les solutions en prenant :

Page 82: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 82/85

Modélisation

FN CRES

82

⎪⎪

⎪⎪

−−=

−=

=

−−−

−−

L

n,2-nn1-n,2-n1n2n2n

n,1-nn1n1n

nn

c xc xdx

c xdx

dx

 

Remarque : Il existe les mêmes problèmes que pour la méthode de Gauss-Jordan quand les

 pivots sont nuls ou proches de zéro.

7.2.2.1   Nombre d'opérations

Pour la méthode de Gauss-Jordan, il faut environ2

n3

opérations. Pour la méthode de Gauss, il

faut environ3

n3

opérations. Ce qui est peu, comparé aux méthodes nécessitant le calcul de

déterminants (méthode de Cramer). Ces approches peuvent donc être utilisées pour des

valeurs de n grandes.

Méthode de Cramer : Si D est le déterminant de la matrice A et si iD est le

déterminant de la matrice obtenue en remplaçant la ième colonne de A par le vecteur B, alors la

solution du système est le vecteur X tel queD

Dx i

i =  

Cette méthode nécessite le calcul de (n+1) déterminants obtenus chacun par (n-1) n!

multiplications, soit (n²-1) n! opérations.

7.3. Méthode du double balayage pour les matrices tridiagonales(Cholesky).

Ces matrices sont par définition relativement creuses. Dans le cas d'un ordre élevé, il y a donc

 beaucoup de zéros. D'un point de vue informatique, on a alors intérêt à stocker une série de

vecteurs représentant chacun une diagonale plutôt qu'une matrice (n × n).

Pour résoudre ce type de matrice, stockée sous forme de vecteurs, on utilise plutôt la méthode

de Gauss que celle de Gauss-Jordan, cette dernière faisant apparaître au cours des calculs desvaleurs non nulles en dehors des trois diagonales.

Soit le système d'équation correspondant à un système matriciel tridiagonal.

⎪⎪

⎪⎪⎪

=+

=++

=++

=+

−−−−−−

nnn1nn

1nn1n1n1n2n1n

2322212

12111

yxaxc

yx bxaxc

yx bxaxc

yx bxa

MMMM  

Page 83: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 83/85

Modélisation

FN CRES

83

qui peut s'écrire

⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢

=

⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢⎢⎢⎢

−−−−

n

1n

2

1

n

1-n

2

1

n

1n

n

1n1n

222

11

y

y

y

y

x

x

x

x

 

a

 b0

0

c

a

0

c0

0

0

0 bac

0 ba

M

M

M

MM

M

LL

MOOOM

OOO

LL

 

 Remarque : 1c et  nb sont nuls.

1er balayage (« aller ») : on crée deux vecteurs A et B . Calcul des termes ii B et A :

0B A 00 ==  

 pour i = 1 à n

i1ii

1iiii

i1ii

ii

aAc

BcyB 

aAc bA 

+

−=

+−=

− 

2ème balayage (« retour ») : on en déduit les valeurs de x

nn Bx =  

 pour i = (n-1) à 1

i1iii BxAx += +  

7.4. Méthodes itératives

Elles sont adaptées aux systèmes linéaires de grande dimension et de matrice peu dense.

7.4.1. PrincipeSoit le système matriciel A X = B. On décompose la matrice A en deux matrices M et N telles

que A = M + N, d'où :

(M + N) X = BM X = -N X + B

X = M-1 (-N X + B)

En partant d'un vecteur X(0) initial, on calcule successivement les nouveaux vecteurs X(1) puis

X(2) …etc par la formule de récurrence :

X(k) = M-1 (-N X(k-1) + B) (attn: en pratique, pas besoin de calculer explicitement 

 M -1

!!) 

On démontre alors que si il y a convergence, on converge sur X, solution du système A X = B

Cette convergence est vérifiée si le rayon spectral de la matrice M-1

  N est inférieur à 1 (lerayon spectral étant le module maximal des valeurs propres = imax λ ).

Page 84: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 84/85

Modélisation

FN CRES

84

 

  Remarques : la condition de convergence n'est pas toujours facile à vérifier. On la

remplacera alors par une condition suffisante telle que ∑≠=

i j j

ijii aa (matrice A à diagonale

dominante)

Le choix des matrices M et N va dépendre de la méthode. Mais M sera choisie de façon à être

facilement "inversible".

On s'arrête quand l'écart des valeurs du vecteur X entre deux itérations devient petit.

7.4.2. Méthode de JacobiDans cette méthode on pose A = D + G + S avec

D matrice diagonaleG matrice triangulaire inférieure

S matrice triangulaire supérieure

S

G

D

 On prend M=D et N=G+S, d'où le schéma itératif :

D X(k) = B – (G + S) X(k-1) 

qui s'écrit sous forme matricielle

)1k (n

2

1

2n1n

n221

n112

n

2

1

)k (n

2

1

nn

22

11

x

x

x

 

0aa

a0a

aa0

 b

 b

 b

x

x

x

 

a00

0

a0

00a

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

M

L

OML

L

MM

L

OMM

L

 

Donc on obtient la formule :

⎥⎥⎥

⎢⎢⎢

⎡−= ∑

≠=

n

i j1 j

 jiji

ii

i 1)-(k )k (xa b

a

1x

(cette formule revient à calculer chaque i kx )( en l’isolant dans la ième

équation où les autres

inconnues sont remplacées par les valeurs calculées à l’itération précédente x j k )1( − ).

 

Il faut pour cela que iia soit non nul quelque soit i (ou que D soit réversible), ce qui peut

généralement être obtenu facilement par inversion des lignes ou des colonnes de la matrice

A).

7.4.3. Méthode de Gauss -SeidelAvec les mêmes notations que précédemment (A = G + D + S), on prend cette fois M= G+D

et N=S.

Le schéma itératif s’écrit donc :

(G + D) X(k) = B – S X(k-1) 

Page 85: cours de modélisation

5/11/2018 cours de mod lisation - slidepdf.com

http://slidepdf.com/reader/full/cours-de-modelisation 85/85

Modélisation

FN CRES

85

soit sous forme matricielle

)1k (n

2

1

n2

n112

n

2

1

)k (n

2

1

nn)1n(n1n

2221

11

x

x

x

 

000

a00

aa0

 b

 b

 b

x

x

x

 

aaa0

aa

00a

−−⎥⎥

⎥⎥

⎢⎢

⎢⎢

⎥⎥

⎥⎥

⎢⎢

⎢⎢

⎥⎥

⎥⎥

⎢⎢

⎢⎢

=

⎥⎥

⎥⎥

⎢⎢

⎢⎢

⎥⎥⎥

⎥⎥

⎢⎢⎢

⎢⎢

MLOM

L

L

MML

OM

M

L

 

Donc on obtient la formule : ⎥⎦

⎤⎢⎣

⎡−−= ∑∑

=+=

1i

1 j

 jij

n

1i j

 jiji

ii

i (k)1)-(k )k (xaxa b

a

1x  

C'est le même schéma que pour la méthode de Jacobi, sauf qu'il faut intégrer les valeurs de

X(k) calculées au fur et à mesure (on calcule chaque i kx )(  par la ième

équation en remplaçant 

les autres inconnues par leur dernière valeur calculée) 

La condition de convergence est la même que pour la méthode de Jacobi. Cette méthode est plus rapidement convergente que la méthode de Jacobi.

7.4.4. Facteur de relaxationIl permet d'accélérer les processus itératifs ou de les rendre convergents (si ils divergent).

Après chaque calcul de i kx )( (par Jacobi ou Gauss-Seidel), on ajuste cette valeur de la

manière suivante : )xx(xx)1k ()k ()1k ()k ( iiii −−

−λ+=  

λ est le facteur de relaxation, = 1 pour une itération normale

< 1 pour une sous-relaxation

> 1 pour une sur-relaxation.

En général :

0 < λ < 1 permet de faire converger un processus divergent

1 < λ < 2 accélérer un processus convergent

λ > 2 divergence

Exemple : la relaxation dans la méthode de Gauss-Seidel donne

⎥⎥⎦

⎢⎢⎣

−⎟⎟ ⎠

 ⎞

⎜⎜⎝ 

⎛ 

−−λ+= −− ∑∑

=+=)1k ((k)1)-(k )1k ()k ( i

1i

1 j

 jij

n

1i j

 jiji

ii

ii xxaxa ba

1

x*x