Cours Edp
Cours Edp
Cours Edp
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
2 Les quations aux drives partielles linaires classiques e e e e 2.1 Equation de transport . . . . . . . . . . . . . . . . . . . . 2.1.1 Modlisation . . . . . . . . . . . . . . . . . . . . . e 2.1.2 Rsolution pour d = 1 . . . . . . . . . . . . . . . . e 2.1.3 Rsolution du cas gnral . . . . . . . . . . . . . . e e e 2.2 Equation de Laplace, quation de Poisson . . . . . . . . . e 2.2.1 Modlisation : membrane lastique . . . . . . . . . e e 2.2.2 Rsultats dunicit . . . . . . . . . . . . . . . . . . e e 2.2.3 Formules de reprsentation . . . . . . . . . . . . . . e 2.2.4 Rgularit des fonctions harmoniques . . . . . . . . e e 2.2.5 Principe du maximum . . . . . . . . . . . . . . . . 2.2.6 Fonctions de Green . . . . . . . . . . . . . . . . . . 2.2.7 Principe de Dirichlet . . . . . . . . . . . . . . . . . 2.3 Equation de la chaleur . . . . . . . . . . . . . . . . . . . . 2.3.1 Modlisation . . . . . . . . . . . . . . . . . . . . . e 2.3.1.1 Equations de raction-diusion . . . . . . . e 2.3.1.2 Equations dadvection-diusion . . . . . . 2.3.2 Calcul dune solution . . . . . . . . . . . . . . . . . 2.3.3 Principe du maximum et unicit . . . . . . . . . . . e 2.3.4 Rgularit sur un domaine born . . . . . . . . . . e e e 2.4 Equation des ondes . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Modlisation . . . . . . . . . . . . . . . . . . . . . e 2.4.1.1 Membrane vibrante . . . . . . . . . . . . . 2.4.1.2 Ondes acoustiques . . . . . . . . . . . . . . 2.4.2 Equation des ondes dans R . . . . . . . . . . . . . . 2.4.3 Mthode des moyennes sphriques (d 2) . . . . . e e 3 2.4.4 Equation des ondes dans R . . . . . . . . . . . . . 2.4.5 Equation des ondes dans R2 . . . . . . . . . . . . . 2.4.6 Mthode de Duhamel pour lquation des ondes non e e 2.5 Classication des e.d.p. linaires dordre 2 . . . . . . . . . e
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . homog`ne e . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3 Equations aux drives partielles dordre un e e 3.1 Introduction et notations . . . . . . . . . . . . 3.2 Etude du cas gnral . . . . . . . . . . . . . . e e 3.3 Cas linaire et quasilinaire . . . . . . . . . . e e 3.4 Lois de conservation scalaires . . . . . . . . . 3.4.1 Modlisation . . . . . . . . . . . . . . e 3.4.2 Etude de lquation de Burgers . . . . e 3.4.3 Solutions faibles et chocs . . . . . . . . 3.4.4 Application : mod`le de trac routier . e
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
47 47 48 50 51 51 52 54 56
4 Introduction aux lments nis ee 59 4.1 Espaces de Sobolev . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2 Formulation faible dun probl`me elliptique linaire . . . . . . . . . . . . . 61 e e e 4.3 Elments nis rectangulaires . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5 Schmas aux dirences nies e e 5.1 Introduction . . . . . . . . . . . . . . . . . . . 5.2 Convergence, consistance et stabilit . . . . . e 5.2.1 Dnitions . . . . . . . . . . . . . . . . e 5.2.2 Condition de Courant-Friedrichs-Lewy 5.3 Analyse de stabilit de von Neumann . . . . . e 5.4 Equations paraboliques . . . . . . . . . . . . . 5.5 Applications . . . . . . . . . . . . . . . . . . . 5.5.1 Equation dadvection-diusion-raction e 5.5.2 Equation de la chaleur . . . . . . . . . 5.5.3 Syst`me de raction-diusion . . . . . e e A Calcul direntiel et intgration dans e e A.1 Notations . . . . . . . . . . . . . . . A.2 Calcul direntiel . . . . . . . . . . . e A.3 Intgration dans Rd . . . . . . . . . . e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 67 69 69 73 74 83 85 85 88 94
Rd 99 . . . . . . . . . . . . . . . . . . . . . 99 . . . . . . . . . . . . . . . . . . . . . 100 . . . . . . . . . . . . . . . . . . . . . 101
B Gomtrie direntielle lmentaire. Intgrale de surface e e e ee e 103 d B.1 Varits dans R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 ee B.2 Espace tangent. Orientation dune varit . . . . . . . . . . . . . . . . . . . 106 ee B.3 Intgrale de surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 e e C Elments danalyse vectorielle 115
Avertissement : ces notes sont un support et complment du cours magistral, e des travaux dirigs et pratiques. Leur contenu nest pas quivalent au cours e e enseign, en particulier les examens et contrles se rf`rent au cours enseign e o ee e uniquement. Bibliographie. Les livres suivants peuvent complter les rappels prsents en annexe : e e e H. Cartan, Cours de calcul direntiel, Hermann 1997 ; e D. Revuz, Mesure et intgration, Hermann 1997. e Des rfrences gnrales sur les quations aux drives partielles sont : ee e e e e e F. John, Partial Dierential Equations, Springer, 4`me dition 1982 ; e e W. Strauss, Partial Dierential Equations : An Introduction, Wiley 1992.
Pour ltude des espaces de Sobolev et la mthode des lments nis on pourra e e ee consulter : H. Brzis, Analyse fonctionnelle, Masson 1987 ; e P.A. Raviart et J.M. Thomas, Introduction a lanalyse numrique des ` e quations aux drives partielles, Masson 1992. e e e Les rfrences suivantes ont t utiles pour la rdaction de ce cours mais leur ee ee e contenu dpasse souvent le cadre dun cours de ma e trise : L.C. Evans, Partial Dierential Equations, AMS 1998 ; G. Folland, Introduction to Partial Dierential Equations, Princeton University Press 1976 ; R. LeVeque, Numerical Methods for Conservation Laws, Birkhuser 1992 ; a I.G. Petrovsky, Lectures on Partial Dierential Equations, Wiley 1954 ; J.C. Strikwerda, Finite Dierence Schemes and Partial Dierential Equations, Wadsworth & Brooks 1989 ; M. Taylor, Partial Dierential Equations, Vol. I, Springer 1996. Un receuil tr`s complet de mod`les mathmatiques en biologie est le livre : e e e J.D. Murray, Mathematical Biology, Springer 1993.
Enn pour les applications en traitement des images on pourra consulter : F. Guichard et J.M. Morel, Image Analysis and Partial Dierential Equations, a para chez Birkhuser. ` tre a Cette liste ne peut tre exhaustive : de nombreux livres traitent des quations e e aux drives partielles et de leurs applications. e e
Chapitre 1 Introduction
1.1 Dnitions e
Soit u une fonction dnie sur Rd , a valeurs dans R et susamment rguli`re pour que e ` e e les expressions qui suivent aient un sens. Une quations aux drives partielles (e.d.p.) pour la fonction u est une relation entre u, e e e les variables x1 , . . . , xd et un nombre ni de drives partielles de u, e e F x1 , . . . , xd , u, D1 u, . . . , Dd u, D1 D1 u, D1 D2 u, . . . , D u, . . . = 0,
o` = (1 , . . . , d ) Nd (voir annexe A.1 pour les notations). u On dit que u est solution de lquation aux drives partielles dans R d si, apr`s e e e e substitution de u et de ses drives partielles, F sannule pour tout (x1 , . . . , xd ) . e e Lordre m N dune quation aux drives partielles est celui de la drive partielle e e e e e dordre le plus lev. e e Un syst`me dquations aux drives partielles est form de plusieurs quations aux e e e e e e drives partielles impliquant une ou plusieurs fonctions inconnues ui . e e Une quation aux drives partielles est linaire si F est linaire par rapport a u et ses e e e e e ` drives partielles. Si m est lordre de lquation aux drives partielles, lquation est de e e e e e e la forme A (x)D u(x) = B(x) ,
||m
1. si u1 et u2 sont deux solutions dune quation aux drives partielles linaire homoe e e e g`ne, alors pour 1 et 2 des rels quelconques, 1 u1 + 2 u2 est aussi solution ; e e 2. si uh est solution de lquation linaire homog`ne et up est solution de lquation e e e e linaire non homog`ne, alors uh + up est solution de lquation compl`te. e e e e
Une quation aux drives partielles dordre m est quasilinaire si F est linaire en toutes e e e e e les drives partielles dordre le plus lev, i.e. dordre m, lquation est alors de la forme e e e e e A (x, u, D u)D u(x) = B(x, u, D u) , avec Nd et || < m .
||=m
Une quation aux drives partielles est semilinaire si lquation est de la forme e e e e e A (x)D u(x) = B(x, u, D u) , avec Nd et || < m .
||=m
La solution dune quation aux drives partielles dordre m dpend en gnral de m e e e e e e fonctions arbitraires de d 1 variables. La solution gnrale dune quation aux drives partielles est celle qui permet de trouver e e e e e toutes les solutions de lquation (sauf des cas de solutions singuli`res) en donnant des e e valeurs particuli`res aux fonctions arbitraires. e Exemples : (i) ux (x, y) = 0 (ii) uxx = 0 (iii) uxx + u = 0 . Pour trouver des solutions particuli`res dune quation aux drives partielles, a partir de e e e e ` la solution gnrale, on va imposer des conditions restrictives sur lensemble des solutions. e e Les contraintes les plus frquentes sont : e 1. conditions initiales : si u est fonction de (x, t) Rd R on donne u(x, t0 ) = 0 (x) ou Dp u(x, t0 ) = p (x) , on parle aussi de conditions de Cauchy ; 2 2. conditions au bord : si u est fonction de x Rd on a trois types de contraintes : conditions de Dirichlet o` u est x sur le bord de : u| = g ; u e conditions de Neumann o` la drive normale de u est x : u e e e conditions de Robin ou mixtes : c(x)u + c(x) du dn =g;
du = g sur ; dn
si g = 0 on a des conditions homog`nes au bord ; e e 3. conditions a linni : si nest pas born on a des conditions de la forme ` u(x) (x) quand |x| ou u 2 < ; 4. conditions sur les interfaces : si = 1 2 , avec 1 2 = 1 2 , et si lon a dtermin u sur 1 et 2 , alors pour pouvoir dnir u sur on a des conditions e e e du , sur 1 2 . sur u, resp. dn Remarques : Les contraintes sont en gnral imposes par la nature du probl`me que e e e e lon essaye de modliser, lquation aux drives partielles et ses conditions restrictives e e e e seront donc a priori cohrentes. e De faon gnrale une quation aux drives partielles ne donne lieu a un probl`me raic e e e e e ` e sonnable que si on lassocie a un certain type de conditions restrictives, par exemple des `
1.1. DEFINITIONS
conditions initiales pour des probl`mes dvolution (quation de la chaleur, quation des e e e e ondes) ou des conditions au bord pour lquation de Laplace. e Exemples : ut (x, t) = 0 sur R u(x, 0) = (x) sur R
2
dans R R u1 (x, y) = 0 + u1 (x, 0) = 0 (x) pour x R Si par ailleurs u1 est solution du probl`me e D2 u1 (x, 0) = 1 (x) pour x R 1 e alors u2 (x, y) = u1 (x, y) + k+1 sin(nx) sinh(ny) est une solution du probl`me de Cauchy n de conditions initiales u2 (x, 0) = 0 (x) et D2 u2 (x, 0) = 1 (x) + nk sin(nx). Les condition initiales sont donc arbitrairement proches mais les solutions u1 et u2 sont arbitrairement loignes, pour y > 0 quelconque (petit). e e
x R,
1.2
Notions de solution
Soit u la solution dun probl`me bien pos. e e Dans ce qui prc`de on na pas prcis la rgularit de u, on a seulement suppos que e e e e e e e u tait susamment direntiable pour que les quations aient un sens, par exemple e e e u Cm () si u est solution dune quation aux drives partielles dordre m . e e e On dit alors que u est une solution au sens fort de lquation aux drives partielles. e e e Dans beaucoup de cas un probl`me nadmet pas de solutions rguli`res : lquation des e e e e ondes avec conditions initiales non continues, les lois de conservation scalaires qui gn`rent e e des ondes de choc. Pour que le probl`me soit bien pos il faut largir lensemble des u, par exemple admettre e e e des solutions non continues. On doit transformer le probl`me de faon a lui garder un sens e c ` pour des fonctions non direntiables, on parle de formulation faible du probl`me. e e Pour une solution forte il y a quivalence entre la formulation faible et le probl`me initial. e e Une solution du probl`me faible est appele solution faible ou solution gnralise. e e e e e Il est souvent plus facile de trouver une solution faible et de montrer a posteriori quelle est rguli`re, cest-`-dire que lon a bien une solution forte. e e a Le choix de solution gnralise (C0 , L1 , . . .) dpend fortement du probl`me trait, de e e e e e e faon gnrale on peut travailler dans le cadre des distributions de L. Schwartz. c e e On verra plus tard dautres espaces de solutions gnralises. e e e d Soit un ouvert de R on dnit e D() = C () = c fonctions de classe C a support compact K , `
cest lensemble des fonctions tests. Lespace des distributions sur est le dual topologique D () des formes linaires continues sur D(). e Exemples : 2. Soit f Lp () = f Lp (K) pour tout compact K , loc alors f D () et, pour tout D() : < f, >=
f (x) (x) dx .
Sans donner de dtails sur la topologie des espaces D() et D (), on rappelle quelques e rsultats de la thorie des distributions : e e - soit (Tk )k une suite de D (), alors Tk T dans D () si pour tout D() : < Tk , > < T, > ; - soit T D (), Nd , on dnit D T par < D T, >= (1)|| < T, D > pour e tout D() ;
- soit S (Rd ) lespace des distributions tempres, i.e. le dual de lespace de Schwartz ee d S(R ) des fonctions a dcroissance rapide : ` e S(Rd ) = f C (Rd ) / , Nd : |x D f (x)| est born . e Pour f S(Rd ) on a f S(Rd ) avec f () = 1 (2)d/2 f (x) ei x. dx .
Rd
Pour T S (Rd ) on dnit sa transforme de Fourier par < T , >=< T, > e e d pour tout S(R ) . En particulier D T () = i|| T () . Exemples : = Application : Considrons loprateur direntiel linaire L = e e e e
||m
e e Si u et v sont dans C () on a une gnralisation de la formule de Green (cf. annexe C) : v L(u) dx = L(v) u dx + M (u, v, n) dS ,
o` L(v) = u
||m
Si L = on a la seconde identit de Green et alors L = . e Si v est a support compact dans on a M = 0 sur , on dit que L est ladjoint de L . ` Pour une distribution T D () on dnit L(T ) par < L(T ), >=< T, L() >, pour e tout D(). On dit que T est une solution fondamentale de ple x0 pour L si L(u) = x0 . o Si a une solution fondamentale, T , on ajoute une solution forte, w Cm (), de lquation ` e homog`ne, L(w) = 0, on obtient encore une solution fondamentale. e Exemple : La fonction u(x1 , x2 ) = 1 pour x1 > x0 , x2 > x0 1 2 0 sinon 2 dans R2 . x1 x2
1.3
1.3.1
Sans rentrer dans les dtails et en ngligeant les constantes dimensionnelles, citons quelques e e quations aux drives partielles cl`bres de la physique : e e e ee u = 0 u f = 0 ut u = 0 utt u = 0 utt uxx + ut + u = 0 quation de Laplace e quation de Poisson e quation de la chaleur, diusion homog`ne e e quation des ondes e quation des tlgraphistes e ee quation de Korteweg-de Vries pour e des vagues sur de leau peu profonde quation de Schr dinger e o
ut + cuux + uxxx = 0 it + = V
et des syst`mes dquations aux drives partielles : e e e e Et = rot(H) quations de Maxwell pour e Ht = rot(E) le champ lectrique E et le champ magntique H e e div(E) = div(H) = 0 Ut + (U. )U U = p div(U ) = 0 quations de Navier-Stokes e dun uide incompressible et visqueux
1.3.2
Dans les mod`les en biologie ou en dynamique des populations on sintresse a la concene e ` tration dune substance chimique ou a la densit dune population (particules, cellules) et ` e son volution au cours du temps. Linconnue u est en gnral fonction de (x, t) Rd R+ . e e e Lquation de Fisher (1937), pour u(x, t) : R R+ R , est un mod`le de dispersion e e en une dimension spatiale dun g`ne favorable dans une population e ut = ku + ru 1 u C ,
le terme de croissance logistique dpend de la constante de reproduction linaire r, de e e la capacit de lenvironnement C, et du coecient de dispersion de la population, k . e Les quations de FitzHugh-Nagumo e ut = uxx + f (u) v vt = vtt + u v pour , , dans R+
utilises pour modliser la transmission dimpulsion nerveuses le long daxones ou des e e ractions chimiques cycliques du type Belousov-Zhabotinsky. e
Les deux quations prcdentes sont des quations de type raction-diusion . e e e e e Les quations de raction-diusion ont t proposes par A. Turing (1952) pour la moe e ee e dlisation de phnom`nes de morphogen`se (i.e. dveloppement des formes/structures). e e e e e Un mod`le dinteraction desp`ces ou de substances chimiques est donn par le syst`me e e e e dquations aux drives partielles e e e ui = div Di t ui + Q i , pour 1 i n ,
o` ui (x, t) reprsente la densit, resp. la concentration, de la substance i . Les matrices de u e e diusion Di et les termes de raction Qi peuvent dpendre de (x, t) et des concentrations e e ui , de faon non linaire. c e Le terme de raction Qi modlise linteraction des substances (inhibition, catalyse) et e e le terme div(Di u) reprsente la diusion de la substance a travers le syst`me. e ` e En fonction du choix de Di et Qi , les concentrations ui peuvent donner lieu a des ` motifs locaux : on modlise ainsi la pigmentation des coquillages, le pelage des animaux e (z`bre, gupard, . . .), cf. Meinhardt (page 94), et des ractions chimiques cycliques, cf. e e e Belousov-Zhabotinsky. En synth`se dimages des chercheurs ont utilis ce mod`le pour crer des textures nae e e e turelles. Lquation de von Foerster e a > 0,t > 0 ua + ut = d(a) u u(a, 0) = f (a) a>0 + u(0, t) = 0 n(a) u(a, t) da t > 0
o` u(a, t) est la densit de population dge a a linstant t , n(a) est le taux de naissances, u e a ` d(a) le taux de dc`s et f (a) la distribution initiale de la population. e e Les graphes de n(a) et d(a) ont lallure indique en dessous. e
n(a) d(a)
PSfrag replacements a a
Exercice : interprtez lquation aux drives partielles, la condition initiale et la condie e e e tion au bord du mod`le. e
1.3.3
Une image numrique scalaire est une matrice I(c, l)0c<N C ,0l<N L o` a chaque ligne l et e u` colonne c on fait correspondre un niveau de gris I(c, l) {0, . . . , 255} , (c, l) est un pixel (angl. picture element) Suite a des probl`mes de capteurs, de transmission ou autres, limage peut tre endomma` e e ge, cest-`-dire que les valeurs en certains pixels ont t modis ou dtruits. Le but des e a ee e e techniques de dbruitage que nous allons exposer est dliminer ce bruit tout en prservant e e e le maximum dinformation contenu dans limage originale. On va modliser une image par une fonction u dnie sur un rectangle ouvert de R2 et a e e ` valeurs dans R, u(x, y) R si lon se restreint a des images en niveaux de gris. ` En traitement du signal on utilise les ltres linaires passe-bas pour enlever le bruit (i.e. e les hautes frquences) dun signal. Le ltre gaussien est un ltre linaire qui est souvent e e utilis car il a une bonne localisation en espace et en frquence. e e 2 1 |x|2 e 2 . Soit uo limage originale bruite et, pour x R2 , posons G (x) = e 2 2 On calcule le rsultat du ltrage de uo par le ltre gaussien : e (G uo )(x) =
R2
G (x y)uo (y) dy .
La valeur de > 0 indique la taille spatiale des structures limines par le ltrage : plus e e est grand, plus on lisse et plus on perd de dtails. e Or on montre que la convolution avec la Gaussienne G consiste a faire voluer limage ` e originale suivant lquation de la chaleur e ut = u u(x, 0) = uo (x, y) . Les proprits de rgularisation de la convolution par G sinterpr`tent donc en termes ee e e de diusion isotrope et le temps t est li a la lchelle spatiale par t = 2 /2 . e` e La diusion isotrope sobtient en prenant une moyenne glissante pondre qui a tendance ee a liminer les pixels bruits mais qui rend ou les contours des objets de limage. `e e Sur la gure de la page suivante on a reprsent quelques tapes de la diusion isotrope e e e dune image bruit. e Pour viter la destruction de contours on prf`re appliquer une diusion anisotrope, ceste ee a-dire que lon va lisser limage partout, sauf au-del` des bords dobjets. Cette mthode ` a e ncessite donc deux tapes. On doit dabord dtecter les contours et autres structures e e e nes de limage, ensuite on va lisser limage en ne considrant que des voisinages de pixels e qui ne traversent pas les bords. On modlise ce comportement par exemple par les quations aux drives partielles non e e e e linaires (1.1) et (1.2), prsentes plus loin. e e e
replacements
PSfrag replacements
Itrations du Laplacien sur une image bruite e e De haut en bas et de gauche a droite : t = i t pour i = 0, . . . , 5 . `
10
En haut a gauche une image bruit u et en haut a droite une dtection de bord par gra` e ` e dient discret | u| . Ici un gradient fort est reprsent en noir et un gradient nul est blanc e e (i.e. inversion vido). e En bas a gauche on a appliqu lquation de la chaleur a u : on calcule la solution de ` e e ` lquation aux drives partielles a linstant t , ceci revient a faire la convolution de u avec e e e ` ` . G 2t En appliquant le gradient a cette image lisse on obtient limage de droite | G2t u| . ` e On remarquera lpaisseur des bords obtenus. e
11
En haut a gauche on a reprsent un lissage obtenu par une technique non linaire mod` e e e e lise par lquation aux drives partielles, dite de la mean curvature motion e e e e u = | u| div t u | u| avec u(x, 0) = uo (x) . (1.1)
En appliquant le gradient on remarque que les bords sont bien localiss : il reste peu de e bruit mais les coins des contours sont arrondis. En bas a gauche un rsultat pour lquation aux drives partielles ` e e e e u = div t u | u| avec u(x, 0) = uo (x) . (1.2)
Limage du gradient montre que le bruit dans les zones homog`nes a compl`tement dise ` e paru mais les bords sont irrguliers. e
12
Interprtation des quations aux drives partielles prcdentes : e e e e e e Dans le rep`re canonique (e1 , e2 ) de R2 on note e Hu (x, y) = uxx (x, y) uxy (x, y) uxy (x, y) uyy (x, y)
la matrice associe a la forme bilinaire symtrique d2 u(x,y) , cest la matrice Hessienne . e ` e e On rappelle que u(x, y) = trace(Hu (x, y)) et u(x, y) = ux (x, y) . uy (x, y)
du = d
u.
et
u = d2 u(, ) = t Hu .
En particulier pour = e1 et = e2 on retrouve les drives partielles dordre un et deux e e habituelles. Supposons que u C (R2 ) et soit c = u(x0 , y0 ), pour (x0 , y0 ) R2 , tel que u1 (c) ne contient pas de points critiques (par le thor`me de Sard ceci est vrai pour presque tout e e c R). En particulier | u(x0 , y0 )| = 0 , le thor`me dinversion locale entra alors que e e ne par le point (x0 , y0 ) il passe une ligne de niveau unique Ic = {(x, y) / u(x, y) = c} . Lensemble dintensit lumineuse constante Ic est appel isophote. Sur la gure suivante e e on a reprsent une isophote ferm. e e e
u PSfrag replacements u<c u (x0 , y0 ) u>c
uy est le vecteur perpendiculaire direct a u en (x, y) . ` ux Dans le rep`re local (, ) on a u = u. = 0 car u est constante sur Ic , tandis que e indique la direction de plus grande croissance de u, u = | u| . Un calcul montre quau point (x, y) Ic : div u | u| = 1 | u| u d2 u u u , | u| | u| u u , | u| | u| .
13
En particulier on remarque que u = d2 u(, ) + d2 u(, ) , on retrouve linvariance du Laplacien et lquation de la chaleur scrit e e ut = u + u . Dans les coordonnes locales (, ) lquation (1.1) scrit e e e ut = u ; tandis que (1.2) devient ut = u . u
On a donc une diusion dans la direction perpendiculaire au gradient, on ne lisse pas les contours dans une image ! Un mod`le qui permet de prendre en compte en chaque point les deux directions et a e t propos par Malik et Perona : ee e u = div g(| u|) t u avec u(x, 0) = uo (x) , (1.3)
o` le coecient de diusion g(r), r 0, est une fonction rguli`re positive dcroissante u e e e avec g(0) = 1 . Dans le rep`re local (, ) lquation (1.3) devient : e e ut = (g(u ) + u g (u )) u + g(u ) u = G (u ) u + g(u ) u , o` on note G(r) = rg(r) lintensit du ux. u e Au points o` | u| = u = 0 on obtient lquation de la chaleur isotrope, pour les autres u e points on obtient une diusion anisotrope pondre par les valeurs de g et G . ee Si G (u ) = 0 on a uniquement une diusion dans la direction perpendiculaire au gradient ; si G (u ) > 0 on ajoute une diusion dans la direction du gradient ; pour G (u ) < 0 on a une diusion inverse instable dans la direction (voir page 25). e 1 Un exemple classique de coecient de diusion est g(r) = 2 , R , + 2 +r 2 r 2 dans ce cas G (r) = 2 et est le seuil qui dcide du signe de G . e ( + r 2 )2 e e e e Note : Les quations aux drives partielles (1.1),(1.2) et (1.3) sont non linaires et leur tude dpasse le cadre de ce cours. Il faut en particulier dnir des solutions faibles e e e adquates pour pouvoir esprer obtenir des probl`mes bien poss. e e e e
14
Equation de transport
Modlisation e
Dans Rd on consid`re un liquide qui volue a la vitesse V (x, t) et dans lequel des particules e e ` de polluant sont introduit. On note u(x, t) la concentration des particules de polluant dans le liquide et f (x, t) la source de polluant. Le champ de vecteurs U (x, t) = u(x, t)V (x, t) reprsente le ux de particules de polluant. e On suppose que les particules de polluant sont simplement entra es par le liquide, i.e. il n ny a pas de diusion du polluant, et u, f et V sont rguliers. e Posons m(t) =
ut (x, t) dx .
Or la variation de m est due a la perte de polluant a travers et a lapparition de ` ` ` particules dans : m (t) =
f (x, t) dx
f (x, t) dx
pour tout Rd , on en dduit la loi dquilibre du polluant e e ut + div(U ) = ut + u.V + u div(V ) = f . Si on suppose que V = c Rd est constante et si (x) est la distribution initiale du polluant on obtient lquation de transport ou dadvection (transport horizontal) e ut + c. u = f (x, t) pour x Rd , t > 0 u(x, 0) = (x) pour x Rd (2.1)
16
Note : dans le cas de transport d a une dirence de temprature on parle de convection ; u` e e des uides chauds (faible densit) montent, tandis que des uides froids (forte densit) e e descendent.
2.1.2
Rsolution pour d = 1 e
ut + cux = 0 pour x R, t > 0 , u(x, 0) = (x) pour x R
du c (x, t) = 0 , , lquation aux drives partielles devient e e e 1 dv cest-`-dire u est constante dans la direction v . a On dnit les droites caractristiques dans le plan tx par x = ct + et on vrie que si u e e e est solution, alors u(ct + , t) = ux c + ut = 0 . t Donc u(x, t) ne dpend que de pour tout (x, t) vriant x ct = . e e Si on pose v =
t x ct = PSfrag replacements x ct =
Donc u(x, t) = F () = F (x ct) est la solution gnrale de lquation aux drives e e e e e partielles et, en utilisant la condition initiale, on obtient la solution du probl`me e u(x, t) = (x ct) pour (x, t) R R+ .
Le domaine de dpendance de u(x, t) par rapport aux valeurs initiales est rduit au point e e . Linuence de la donne en est rduite au points de la droite caractristique. e e e Le graphe de u a linstant t est celui de translat de ct . ` e Interprtation : Considrons un tube dans lequel coule de leau a la vitesse c > 0 et qui e e ` contient des traces de polluant. Dapr`s ce qui prc`de la distribution de polluant que lon e e e a en t = 0 se retrouve inchange a linstant t1 > 0, a une translation de +ct1 pr`s . e ` ` e
PSfrag replacements t=0 c x = ct1 t = t1 c
x=0
x=0
17
2.1.3
En sinspirant du cas monodimensionnel on consid`re les droites caractristiques dans e e lespace des (x, t) Rd+1 : c une droite caractristique passant par le point (x, t) et de direction v = e Rd+1 est 1 paramtrise par (x + s c, t + s) , pour s R. e e On pose z(s) = u(x + s c, t + s) et, si u est solution de (2.1), on a z (s) = f (x + s c, t + s) , do` u
0 0
z (s) ds =
t t
f (x + s c, t + s) ds ,
u(x, t) = (x t c) +
f (x + (s t) c, s) ds ,
pour (x, t) Rd R+ .
18
2.2
2u (x) = x2 i
Cest lquation de Poisson qui modlise des phnom`nes dquilibre : potentiel lectroe e e e e e statique, membrane en quilibre, champ gravitationel,. . . e Si f = 0 on a lquation de Laplace et on dit que u est une fonction harmonique . e
2.2.1
On consid`re une membrane lastique qui sidentie a une rgion plane R2 si on e e ` e napplique aucune force. Si on applique une force f , normale a R2 , le membrane se dforme ` e et prend une position dquilibre M . e En supposant des petites dformations on note u(x) M la position dquilibre du point e e x = (x1 , x2 ) , donc M = u() .
u
u(x)
PSfrag replacements x2
x1
1 + | u|2 1 dx , avec R . +
Comme |ux1 | et |ux2 | sont petits, on a 1 + | u|2 1 + | u|2 /2 . Le travail de la force f est donn par le dplacement des points x , on en dduit e e e lexpression de lnergie potentielle de la membrane a lquilibre : e ` e E(u) =
1 | u|2 + f u 2
dx .
Posons () = E(u + ), o` u est une dformation admissible de M , i.e. compatible avec u u e les contraintes imposes sur u . Comme la membrane est a lquilibre, E(u) est minimal, e ` e donc 0 = (0) = u. u + f u dx .
du u d = 0 . dn
()
19
Des contraintes sur le bord de la membrane M donnent, a travers lintgrale curviligne ` e dans (), des contraintes sur u : si le bord de la membrane est x, u(x) = b(x) pour x , on a u| = 0 et () est e vri ; on obtient le probl`me de Dirichlet e e e 1 u = f sur ; u= b sur si le bord de la membrane nest pas x, u est non nul sur et () est vri si u est e e e solution du probl`me de Neumann avec conditions au bord homog`nes e e u = 1 f sur ; du = 0 sur dn
si une force f1 , normale a R2 , agit sur le bord de la membrane, lintgrale curviligne ` e du dans () devient + f1 u d, et on a le probl`me de Neumann e dn u = 1 f sur . du 1 = f1 sur dn
2.2.2
Rsultats dunicit e e
` Soit un ouvert born de Rd , de bord rgulier et u C2 () . A partir des formules e e de Green on obtient, en prenant v = u, resp. v = 1 uu dx =
du dS dn
| u|2 dx .
(2.2)
u dx =
du dS dn
(2.3)
Si u est harmonique dans on obtient grce a (2.2) la proposition suivante. a ` Proposition 2.2.1 e Une solution u C2 () du probl`me de Dirichlet u = f sur est unique. u = g sur
u = f sur 2 e est unique a une ` Une solution u C () du probl`me de Neumann du = g sur dn constante pr`s. e
CHAPITRE 2. LES EQUATIONS AUX DERIVEES PARTIELLES LINEAIRES CLASSIQUES
20
Remarque : En utilisant (2.3) on constate quune solution du probl`me de Neumann e existe seulement si f (x) dx = g(x) dS .
2.2.3
Formules de reprsentation e
Fixons x0 Rd , le Laplacien tant invariant par translations et rotations, on se propose e de chercher une fonction harmonique v telle que
d
1 2
v(x) = (r)
pour x Rd avec
r = |x x0 | =
i=1
(xi x0 )2 i
si d = 2 C ln(r) + C , avec (C, C ) R2 . dont la solution est (r) = C 2d r + C si d > 2 2d Donc v(x) = (|x x0 |) est harmonique pour x = x0 . Soit u C2 () et v(x) = (|x x0 |), o` x0 x, en posant = \ Bd (x0 , ) pour u e > 0, on a du dv u dS . vu dx = v dn dn En faisant tendre vers 0 on montre que vu dx =
d1 (r) = 0 , r
du dv u dn dn
dS + Cd u(x0 ) .
(2.4)
1 , C = 0 et on pose d 1 ln |x x0 | 2
si d = 2
K(x, x0 ) est une solution fondamentale de ple x0 de lquation de Laplace , o e en eet x K = x0 au sens des distributions. On a, grce a (2.4), pour tout x0 , a ` u(x0 ) =
K(x, x0 )u(x) dx
K(x, x0 )
dSx ,
(2.5)
dont on dduit la e
21
x2
x1
ln
x2 + x 2 1 2
dSx .
2.2.4
Soit w C2 () et w harmonique dans , alors G(x, x0 ) = K(x, x0 ) + w(x) est encore une solution fondamentale de ple x0 de lquation de Laplace. Grce a (2.5) on a o e a ` u(x0 ) =
G(x, x0 )u(x) dx
G(x, x0 )
dSx ,
(2.6)
et on en dduit le e Thor`me 2.2.1 (Loi de la moyenne arithmtique de Gauss) e e e Soit u harmonique dans Bd (x0 , r) et continue sur Bd (x0 , r), alors u(x0 ) = 1 d r d1 u(x) dS =
S d1 (x0 ,r)
1 d rd
u(x) dx .
Bd (x0 ,r)
Remarques : 1) On peut montrer rciproquement quune fonction u C2 () qui vrie la proprit de e e ee la moyenne pour toute Bd (x0 , r) est harmonique dans . 2) De faon plus gnrale, si u(x) 0 dans Bd (x0 , r) on a c e e u(x0 ) 1 d r d1 u(x) dS
S d1 (x0 ,r)
22
Thor`me 2.2.2 e e Si u C0 () vrie la proprit de la moyenne pour toute Bd (x0 , r) , alors u C (). e ee Remarques : 1) Noter que u nest pas ncessairement continue sur . e 2) On montre quune fonction harmonique sur est en fait analytique.
2.2.5
Principe du maximum
Thor`me 2.2.3 (Principe faible) e e Si u C2 () C0 () et u 0 sur , ouvert born de Rd , alors e max u = max u .
Thor`me 2.2.4 (Principe fort) e e Si u C2 () C0 () et u = 0 sur , ouvert born connexe de Rd , alors e soit soit u est constante dans , min u < u(x) < max u , pour tout x .
2.2.6
Fonctions de Green
Dfinition 2.2.1 e On dit que G(x, x0 ) est une fonction de Green pour le probl`me de Dirichlet pour e lquation de Laplace dans , si e (i) G(x, x0 ) = K(x, x0 ) + w(x, x0 ) pour x , x0 et x = x0 (ii) G(x, x0 ) = 0 pour x , x0 o` w C2 () et w = 0 dans ; K est dnie page 20. u e Par dnition G est une solution fondamentale et (2.6) donne, pour tout x0 e u(x0 ) =
u(x)
dG (x) dS . dn
(2.7)
u = 0 sur . u = g sur
Pour construire G il faut dterminer w, ce qui revient a rsoudre un nouveau probl`me e ` e e de Dirichlet pour lquation de Laplace dans . e Il est en gnral dicile de trouver une expression explicite pour G. On va prsenter un e e e cas o` la gomtrie de permet de dterminer G. u e e e
23
Fonction de Green pour une boule. Soit = Bd (0, a) et x0 Bd (0, a) x. e On montre que = S d1 (0, a) est le lieu de points x Rd qui vrient e o` x0 = u a |x0 |
2
|x x0 | a = |x x0 | |x0 | a r. |x0 |
On construit la fonction de Green G(x, x0 ), a partir de K(x, x0 ) et K(x, x0 ), en remar` quant que K(x, x0 ) est harmonique dans Bd (0, a) et de classe C2 pour tout x = x0 . Pour d = 2 on consid`re la fonction de Green : e 1 G(x, x0 ) = ln |x x0 | + ln 2 Pour d 3 on a la fonction de Green : 1 G(x, x0 ) = (2 d)d a |x0 |
ln |x x0 |
2
1 |x x0 |d2
a |x0 |
1 |x x0 |d2
Grce a (2.7) on obtient la a ` Proposition 2.2.3 (Formule de Poisson) Soit u C2 (Bd (0, a)) et u = 0 dans Bd (0, a), alors pour tout x0 Bd (0, a) u(x0 ) =
S d1 (0,a)
o` H(x, x0 ) = u
En tudiant les proprits de H on montre rciproquement la e ee e Proposition 2.2.4 Soit g C0 (S d1 (0, a)), alors la fonction u dnie par e g(x0 ) si |x0 | = a u(x0 ) = H(x, x0 ) g(x) dSx si |x0 | < a
S d1 (0,a)
Remarque : Un deuxi`me cas o` la gomtrie de permet de dterminer facilement G e u e e e est celui dun demi-espace = Rd1 R+ . On a = (Rd1 R+ )={x Rd / xd = 0} = Rd1 . Pour x0 Rd1 R+ on obtient les quivalents des propositions 2.2.3 et 2.2.4 avec e u(x0 ) =
Rd1
H(x, x0 ) u(x) dx ,
o` H(x, x0 ) = u
24
2.2.7
Principe de Dirichlet
Soit un ouvert, born de Rd . Pour tout v dans A = v C2 () / v = g sur e on dnit la fonctionnelle dnergie e e 1 | v|2 v f 2
E[v] =
dx ,
o` f C0 () . On a le rsultat de calcul des variations suivant (cf. section modlisation) : u e e Proposition 2.2.5 e Soit u C2 (), on a quivalence entre (i) u est la solution unique de u = f sur u = g sur
25
2.3
Equation de la chaleur
u ut (x, t) ku(x, t) = (x, t) k t
d
i=1
pour u dni sur Rd R et k > 0 . e + Cest lquation de la chaleur qui modlise des phnom`nes dvolution : diusion de chae e e e e leur, rpartition de substances chimiques, mlange desp`ces,. . . e e e Remarques : 1) Un changement dchelle, t = k t , transforme lquation aux drives partielles en e e e e 2 ue = u , il sut donc dtudier le cas k = 1 . Comme k [m /s], on obtient une quation e e t aux drives partielles normalise sans dimensions. e e e 2) Poser t = t change compl`tement lquation aux drives partielles : on ne peut pas e e e e inverser le temps. 1 2 Contre-exemple (Petrovsky) : pour tout n N la fonction u(x, t) = sin(nx) en kt n est une solution de lquation de la chaleur sur R R. e 1 On a u(x, 0) = sin(nx) qui est arbitrairement petit pour n grand. n 1 2 En t = + > 0, la solution u(x, +) = sin(nx) en k est borne par u(x, 0) . e n 1 2 ` Par contre, en t = , on a u(x, ) = sin(nx) e+n k qui explose par rapport a n u(x, 0) . 1 Le probl`me ut = kuxx sur R R, avec u(x, 0) = sin(nx), est instable donc mal pos. e e n Lquation de la chaleur modlise des phnom`nes dvolution irrversibles. e e e e e e 3) Poser (x, t) = (ax, a2 t) prserve lquation aux drives partielles et la quantit |x|2 /t . e e e e e
2.3.1
2.3.1.1
Modlisation e
Equations de raction-diusion e
Nous allons modliser le comportement de diusion dune population (cellules, insectes) e ou de particules (substances chimiques). On suppose lexistence dune source de particules (naissance, resp. dc`s, dinsectes). e e Soit (x, t) R+ , avec ouvert born de Rd et rgulier. On note e e u(x, t) la fonction de densit des particules (la concentration), e e e q(x, t, . . .) le taux de cration net de particules ( naissances moins dc`s ) et e F (x, t, . . .) la densit du ux de particules, cest-`-dire F (x, t).n est le ux de particules e a (par unit de temps) a travers un lment de surface plane, perpendiculaire a n en x et e ` ee ` daire 1 (cf. annexe C).
26
On suppose pour la suite que u et F sont rguliers et on consid`re O de bord rgulier. e e e La variation de masse dans O est due a la cration/destruction de particules dans O et ` e au ux de particules a travers O ` t u(x, t) dx =
O O
q(x, t) dx
Pour exploiter le mod`le on va se donner F et q , deux cas sont prsents : e e e Si u(x, t) est la temprature, alors la quantit de chaleur dans a linstant t est donne e e ` e par simplier, que c = 1 . Loi de Fourier : la chaleur va des rgions chaudes vers les rgions froides a une vitesse proe e ` portionnelle a la variation de temprature. ` e
On suppose de plus que lon ne peut perdre de la chaleur que par , on a donc F = k u , o` k(x, t) est la conductivit de chaleur, et q = 0 . u e Do` ut = div(k u) = u k. u + ku et, si k est constant, ut = k u .
On a obtenu lquation de la chaleur. e Si u(x, t) reprsente la concentration de particules on a la e Loi de Fick : les particules vont des hautes densits vers les faibles densits et leur ux est e e proportionnel a la variation de la densit. ` e On obtient encore F = D u , o` D est le coecient de diusion, et do` lquation u u e de raction-diusion e ut = div(D u) + q . Exemple : Soient A et B des substances chimiques qui ragissent suivant la loi A+B e 2B + R . On note a, resp. b, la concentration de A, resp. B . Si on a un mlange parfait de A et B les concentrations vrient le syst`me direntiel e e e e ordinaire at = ab . bt = ab Si A et B nont pas t mlangs et si on veut modliser leur diusion, on obtient le ee e e e mod`le de raction-diusion e e at = D1 a ab . bt = D2 b + ab
27
2.3.1.2
Equations dadvection-diusion
On consid`re le mouvement de particules (atomes, insectes) sur une grille discr`te monodie e mensionnelle xi = i ( > 0). On suppose que ces particules peuvent changer de position tous les > 0 et quelles bougent dau plus . On note tn = n . On note un la probabilit dtre en xi a linstant tn , avec e e ` i
xi
un i
(0)
u(tn , x) dx = 1 .
R
De plus on note pn la probabilit de passer de la position xi a la position xi+1 entre tn et e ` i n tn+1 , et qi la probabilit de passer de la position xi a la position xi1 entre tn et tn+1 . e `
PSfrag replacements
tn+1 tn
x xi1 xi xi+1
La probabilit pour une particule dtre en xi a linstant tn+1 est donn par lquation de e e ` e e
Chapman-Kolmogorov
n n un+1 = pn un + qi+1 un + (1 pn qi )un , i1 i1 i+1 i i i
o` les trois termes du membre de droite reprsentent respectivement la probabilit daru e e river en xi depuis xi1 , depuis xi+1 ou de rester en xi , entre tn et tn+1 . Si on rcrit lquation sous la forme ee e 1 un+1 un i i = 2
n n pn qi+1 pn qi1 i1 un i+1 un i1 i+1 n n n pn + qi+1 2 n 1 pn + qi1 2 n pn + q i 2 n i1 + 2 ui1 + i+1 ui+1 2 i ui 2 2 2 n pn q i i c(x, t) n pn + q i 2 i d(x, t) 2
et si on suppose que
et
quand et tendent vers 0, on obtient lquation de Fokker-Planck : e u 2 = c(t, x) u + 2 d(x, t) u , t x x o` c(x, t) est la vitesse dadvection (ou convection) et d(x, t) le coecient de diusion. u e Note : les limites ci-dessus se justient dans le cadre de la thorie des processus de Markov. n Si c > 0, i.e. pn > qi , les particules ont un mouvement vers la droite, si c < 0 on a un i mouvement vers la gauche. n Si d est grand, i.e. pn + qi grand ( 1), les particules ont tendance a se dplacer. ` e i
28
n Si pn = qi = 1/2 alors c(x, t) = 0 et d(x, t) = D est constante, on obtient lquation de la e i chaleur ut = D uxx .
Si c(x, t) = C et d(x, t) = D on a une quation dadvection-diusion : ut Duxx +Cux = 0. e e ` Interprtation : considrons un tube dans lequel coule de leau a la vitesse C > 0 et e contenant des traces de polluant. La distribution initiale de polluant est, a linstant t1 > 0, ` transport par advection de Ct1 et modie par diusion. e e
PSfrag replacements t=0 C C x = Ct1 t = t1
x=0
x=0
Remarque : Si on modlise des phnom`nes thermiques on parle de convection et donc e e e dquations de convection-diusion : les mouvements de masses dair dans latmosph`re e e et les courants ocaniques dus a des dirences de temprature inuencent le climat de e ` e e la terre. Des phnom`nes de convection thermique sous lcorce terrestre et dans le noyau e e e sont responsables de la tectonique des plaques et de la cration du champ magntique e e terrestre.
2.3.2
On consid`re le probl`me e e
On va dterminer, de faon formelle, une solution en utilisant la transforme de Fourier e c e 1 u(x, t) ei x. dx . de u par rapport aux variables despace x : v(, t) = (2)d/2 Rd On obtient lquation direntielle ordinaire en v : e e do` u v(, t) = f () e|| u(x, t) = 1 (2)d
2t
et f () e|| t ei x. d =
Rd Rd
|xy|2 1 e 4t . (4t)d/2 2
1 (2)d/2
K(x, y, t) f (y) dy ,
o` K(x, y, t) = u Remarques :
Rd
|x|2 1 e 4t I{t>0} est une solution fondamen(4t)d/2 d tale, de ple (0, 0) R R , de lquation de la chaleur. o e |y|2 1 2. On a K(x, y, t) = G2t (x y) avec G (y) = e 22 , la Gaussienne centre e (2 2 )d/2 dcart type . e
29
PSfrag replacements
x2
x1
K poss`de les proprits suivantes : e ee a) K(x, y, t) est de classe C pour x, y Rd et t > 0 . b) Pour x, y Rd et t > 0 , c)
Rd
x K(x, y, t) = 0 . t
K(x, y, t) dy = 0
|yx|>
Grce a ces proprits on dmontre que u est bien une solution, on a le a ` ee e Thor`me 2.3.1 e e e e Soit f C0 (Rd ), i.e. f continue et borne, alors la fonction dnie par b u(x, 0) = f (x) u(x, t) = 1 (4t)d/2 e
Rd
|xy|2 4t
e est C sur Rd R , vrie ut = u sur Rd R et est continue sur Rd R+ . + + Remarques : 1. Si f est born on a pour tout x Rd et t > 0 : inf f u(x, t) sup f . e
Rd Rd
2. On peut assouplir les hypoth`ses du thor`me : e e e 2 soit f une fonction mesurable et supposons que pour tout y Rd , |f (y)| M ea|y| , avec a, M des constantes xes. Alors la fonction u, dnie dans le thor`me, vrie e e e e e 1 ut = u et est C en tout (x, t) Rd ]0, 4a [ . Si f est continue en x0 on a en plus lim u(x, t) = f (x0 ) .
(x,t)(x0 ,0)
Donc mme avec des donnes initiales non continues la solution devient, a linstant e e ` t > 0, de classe C : lquation de la chaleur a un eet rgularisant. e e
30
3. Lexemple de Tikhonov montre que la solution u du thor`me nest pas unique. e e + (k) 2 (t) 2k e1/t pour t > 0 2 Soit (t) = et, pour (x, t) R , posons u(x, t) = x . 0 pour t 0 (2k)! k=0 Alors u C (R2 ) , vrie ut = uxx sur R R et, pour tout x R, u(x, 0) = 0 . e 4. Le thor`me prcdent permet de considrer la fonction u(x, t) = e e e e e
|x|2 1 e 4kt , (4kt)d/2 d dnie pour (x, t) R R+ , comme solution au sens des distributions du probl`me : e e
sur Rd R + pour x Rd
5. La valeur en (x, t) de u, dnie dans le thor`me prcdent, dpend, pour t > 0, de e e e e e e tous les points y Rd . Rciproquement la valeur de f en y0 Rd aecte a linstant e ` t > 0 (petit) les valeurs de u en tout point x. Les eets de lquation de la chaleur e voyagent a vitesse innie . ` On sintresse donc plutt a la vitesse de diusion dune certaine proportion de e o ` la concentration, resp. chaleur, initiale. Considrons une diusion modlise par e e e x2 1 4kt u(x, t) = , pour (x, t) R R+ . Soit ]0, 1[ donn, a linstant t > 0 e ` e 2 kt la quantit de particules se trouve dans lintervalle [x , x ] et e 1 = 2 kt
x
e
x
x 4kt
1 dx = kt
e
0
x 4kt
2 dx =
x 2 kt
ez dz .
Lintgrale de droite tant constante, on trouve que x est proportionnel a 2 kt . e e ` On dit que t = x2 /k est le temps de diusion. On peut comparer le temps de diusion au temps dadvection t = x/c , voir p. 15 .
2.3.3
On va traiter deux cas, celui dun domaine spatial born et celui de lquation de la chaleur e e d sur R . Cas dun domaine spatial born. e On a = avec = [0, T ] {t = 0} et = {t = T } .
PSfrag replacements T
Rd
31
Thor`me 2.3.2 (Principe du maximum) e e Soit u C0 () tel que ut , uxi xj existent et soient continues dans . Si ut u 0 dans , alors max u = max u .
Thor`me 2.3.3 e e Soit u C0 () tel que ut , uxi xj existent et soient continues dans . Alors u est dtermin de faon unique dans si on se donne ut u sur et u sur . e e c Corollaire 2.3.1 Soit C0 (]0, T ]), f C0 () et g C0 ( [0, T ]). Il existe au plus une solution u du probl`me e ut u = sur ]0, T ] u = f sur {t = 0} , u = g sur [0, T ]
Cas o` le domaine spatial est Rd . u Thor`me 2.3.4 (Principe du maximum) e e Soit f C0 (Rd ) et soit u C0 (Rd [0, T ]) tel que ut , uxi xj existent et soient continues dans Rd ]0, T [ , si de plus (i) ut u 0 (ii) u(x, 0) = f (x)
2
sur Rd [0, T ] .
Corollaire 2.3.2 Soit C0 (Rd [0, T ]), f C0 (Rd ). Il existe au plus une solution u C0 (Rd [0, T ]) du probl`me e ut u = sur Rd ]0, T [ u = f sur Rd {t = 0} qui vrie |u(x, t)| M ea|x| sur Rd [0, T ], pour des constantes M , a > 0 . e Remarques : La solution obtenue dans la section 2.3.2 est unique si f est born et si lon se restreint e aux solutions bornes. e
2
32
Si f est major par une expression du type M ea|x| et si on se restreint aux solutions e vriant le mme type de majorations, on obtient encore une solution unique (cf. remarque e e 3, section 2.3.2). 2 Lexemple de Tikhonov ne peut pas vrier une majoration du type M ea|x| et permet de e construire des solutions tr`s croissantes a partir de nimporte quelle condition initiale f . e `
2.3.4
Soit un ouvert born de Rd avec rgulier, on pose = ]0, T [ . e e On suppose que u, ut et les uxi xj existent et sont continues dans et que ut u = 0 dans . Pour v C2 () quelconque on a 0=
v(ut u) dxdt =
u(vt + v) dxdt
v
0
du dv u dSx dt . dn dn
En particulier on pose v(x, t) = K(x, y, T +t) pour y x et > 0 . Alors v C2 () e et vt + v = 0, donc (2.8) devient K(x, y, )u(x, T ) dx =
T
K(x, y, T + )u(x, 0) dx
+
0
Si on fait tendre vers 0 le membre de gauche de (2.9) tend vers u(y, T ), en eet il est solution de w (y, ) y w(y, ) = 0 avec w(y, 0) = u(y, T ) (mme dmonstration que e e pour le thor`me 2.3.1). e e Si on pose K(x, y, s) = 0 pour s 0, alors K(x, y, s) est C pour x, y Rd , x = y et s R. On obtient u(y, T ) =
+
R
33
2.4
Dans cette section on va tudier lquation aux drives partielles linaire dordre deux e e e e e utt (x, t) c2 u(x, t) = pour u dni sur Rd R et c > 0 . e Cest lquation des ondes qui modlise des phnom`nes dvolution : corde ou membrane e e e e e vibrante, ondes acoustiques, ondes lectromagntiques, ondes sismiques, . . . e e On note u= 2 t2
d
2u (x, t) c2 2 t
i=1
2u (x, t) = 0 , xi 2
i=1
2 xi 2
u le dAlembertien de u.
Remarque : Une inversion et translation du temps, t = t0 t, ne change pas lquation e aux drives partielles, on va restreindre ltude a t 0 . e e e `
2.4.1
2.4.1.1
Modlisation e
Membrane vibrante
Considrons une membrane lastique en vibrations M . Pour x = (x1 , x2 ) R2 et e e t > 0, on note u(x, t) la position de la membrane au point x a linstant t : Mt = u(, t) . ` Alors ut (x, t) est la vitesse du dplacement (vertical) de la membrane en x a linstant t et e ` utt (x, t) est lacclration. ee On suppose des petits dplacements transversaux et on introduit les constantes R , e + la densit de la membrane, et R , le coecient de tension. On obtient lexpression e + suivante pour lnergie potentielle du syst`me a linstant t , voir aussi page 18, e e ` Ep = et lnergie cintique du syst`me est e e e Ec = 2 ut (x, t)2 dx .
| u(x, t)|2 dx
Dapr`s le principe de Hamilton de la mcanique, la membrane passe de ltat en t = 0 e e e a ltat en t = T , T > 0, par un point critique de la fonctionnelle ` e
T
A(u) =
0
(Ec Ep )dt =
1 2
T 0
u2 | u|2 dxdt . t
i.e. la variation premi`re de A est nulle. Posons () = A(u + ), o` la perturbation e u u admissible u est rguli`re et vrie u(x, 0) = u(x, T ) = 0 pour tout x . On a e e e
T
0 = (0) =
0
(ut ut u. u) dxdt ,
34
(utt + u) dxdt u
du u dSx dt = 0 . dn
()
Si lon choisit u tel que u(x, t) = 0 pour tout x et t [0, T ] , on obtient grce au a 2 premier terme de la relation () lquation des ondes dans R : e utt = c2 u , o` c = u . du (x, t) = 0 pour dn
Le second terme de la relation () impose des conditions au bord : si la membrane est xe au bord, u(x, t) = 0 ; si la membrane est libre, e (x, t) R ; + Pour compl`tement modliser le mouvement de la membrane on a des conditions initiales : e e la position u(x, 0) et la vitesse ut (x, 0) , donnes pour x et t = 0 . e 2.4.1.2 Ondes acoustiques
On modlise un uide (liquide ou gaz) par la densit des particules u(x, t) et la vitesse e e 3 V (x, t) pour x R et t > 0 . La densit du ux de particules uV [particules/m2 s] e (cf. annexe C) est ici interprt comme tant la quantit de mouvement ([kg m/s]) des ee e e particules par unit de volume. e
PSfrag replacements x
F (x )
p(x, t) n(x)
On suppose que u et V sont rguliers. Soit un ouvert born quelconque de R3 , de bord e e rgulier ; pour x on note n(x) le vecteur normal unit extrieur. e e e La loi de conservation de masse donne u(x, t) dx = u(x, t)V (x, t).n(x) dSx , t do`, u
ut dx =
div(uV ) dx ,
Ecrivons maintenant la loi de conservation de la quantit de mouvement pour chacune e des trois composantes de uV : t
uVi dx =
uVi V.n(x) dS
pni (x) dS +
uFi dx ;
()
35
o` i {1, 2, 3} , p(x, t) ([N/m2 ] = [kg/ms2 ]) est la pression et F (x, t) la somme des forces u extrieures qui agissent au point x (e.g. gravitation). e Le membre de gauche reprsente la variation de la quantit de mouvement dans , le e e premier terme du membre de droite reprsente le ux de la quantit de mouvement a e e ` travers , le second la pression sur la paroi et le troisi`me lapport des forces externes. e En appliquant le thor`me de la divergence on a, pour 1 i 3 : e e
()
Les quatres quations aux drives partielles (2.11),(2.12) sont compltes par la loi dtat : e e e ee e p(x, t) = p(u(x, t)) ; la pression est fonction croissante de la densit du uide. e Noter que le terme dadvection (V. )V = DV V est non linaire et quil est en partie a e ` lorigine de phnom`nes complexes, et intressants, en mcanique des uides. e e e e e e Approximation acoustique : Pour simplier notre mod`le on va ngliger les forces externes, donc F = 0. De plus, la propagation du son dans un gaz (e.g. lair) est le rsultat de e vibrations que lon va supposer faibles : V est petit et le terme dadvection dans (2.12) est ngligeable. Finalement, on suppose que la densit et la pression scartent peu de e e e quantits nominales u0 et p0 : e u(x, t) = u0 + u(x, t) , u << u0 et p(x, t) = p0 + p(x, t) , p << p0 .
On linarise les quations de la dynamique des uides autour de (u, p, V ) = (u0 , p0 , 0) . e e u dp (u0 ) u et lquation (2.11) devient : e + u0 div(V ) = 0 , Alors p = du t V dp tandis que lquation (2.12) scrit : e e u0 = (u0 ) u . t du En drivant la premi`re quation par rapport a t et en utilisant la seconde, on obtient e e e ` 2u = u0 div 2 t V t 1 dp (u0 ) u0 du u = dp (u0 ) . u du
= u0 div On pose c =
dp 2u (u0 ) et u vrie donc lquation des ondes dans R3 : e e = c2 . u 2 du t La propagation des ondes acoustiques se fait par augmentation et diminution de la densit e de particules. Pour lair la vitesse des ondes est c = 340 m/s .
36
Remarque : Les ondes sur une membrane ou corde vibrante sont dites ondes transversales car elles se propagent dans la direction perpendiculaire au mouvement des particules. Les ondes acoustiques sont des ondes longitudinales car elles se propagent dans la mme e direction que les particules.
Direction de propagation des ondes
PSfrag replacements
Ondes transversales
Ondes longitudinales
2.4.2
On consid`re le probl`me e e utt (x, t) = c2 uxx (x, t) sur R R , c R + + u(x, 0) = f (x) pour x R ut (x, 0) = g(x) pour x R Grce au changement de variables a
(2.13)
= x + ct on obtient lquation u (, ) = 0 . e = x ct
Le domaine de u tant connexe on trouve u(, ) = F () + G(), o` F et G sont des e u 2 fonctions dune variable relle de classe C , et on a e u(x, t) = F (x + ct) + G(x ct). La solution gnrale de lquation des ondes dans R est donc obtenue par superposition e e e de v, solution de vt cvx = 0 : v(x, t) = F (x + ct), et de w, solution de wt + cwx = 0 : w(x, t) = G(x ct).
PSfrag replacements F (x) = G(x) u(x, t) = F (x + ct) + G(x ct) u
x ct +ct
Le graphe de u dans le plan xu montre deux ondes, solutions de deux quations de e transport, qui se propagent, sans changer de forme, vers la gauche pour v et vers la droite pour w . En utilisant les conditions initiales de (2.13) on trouve F (x) = 1 f (x) + 2 2c
x
g(y) dy + et G(x) =
0
f (x) 1 2 2c
x 0
g(y) dy .
37
Si f C2 (R) et g C1 (R), le probl`me (2.13) admet une solution unique de classe C2 sur e R R+ , donne par la formule de dAlembert e u(x, t) = 1 1 f (x + ct) + f (x ct) + 2 2c
x+ct
g(y) dy .
xct
(2.14)
Remarques : 1) La valeur de u(x, t) dpend des valeurs de f et g sur lintervalle [x ct, x + ct], cest e le domaine de dpendance. e Rciproquement, le cne dinuence du point (, 0) est lensemble (x, t) / |x| ct . e o Linformation voyage a vitesse nie c. `
PSfrag replacements t (x, t) t x + ct = x ct =
Domaine de dpendance e
2) La fonction (2.14) est au plus aussi rguli`re que les conditions initiales f et g . e e 3) Pour c > 0 on consid`re lquation des ondes utt = c2 uxx dans R R () . e e + On note (ABCD)/c le paralllogramme form par les points A, B, C et D et dont les e e PSfrag replacements cts se trouvent sur des droites caractristiques de () . oe e
t x + ct = 3 A B C x (1 , 0) (2 , 0) (3 , 0) (4 , 0) D x ct = 2 x + ct = 4 x ct = 1
e e e e Soit u C2 (R R ), alors u vrie lquation aux drives partielles () si et seulement + si pour tout (ABCD)/c R R+ : u(A) + u(C) = u(B) + u(D) .
38
2.4.3
Soit h C0 (Rd ), on dnit la moyenne sphrique de h sur la sph`re S d1 (x, r) par e e e 1 d h(x + r) dS .
||=1
La fonction Mh est ainsi dnie pour tout r R et on a Mh (x, r) = Mh (x, r) . e d Si h C (R ), alors Mh C (Rd+1 ). r 1 En particulier si h C2 (Rd ) on calcule d1 Mh (x, ) d . Mh (x, r) = d1 x r r 0 On en dduit lquation de Darboux, vrie par toute moyenne sphrique dune fonction e e e e e 2 h de classe C d1 2 r 2 + r r Mh (x, r) = x Mh (x, r) (2.15) Mh (x, 0) = h(x) Mh (x, 0) = 0 r
Soit u C2 (Rd R+ ) solution du probl`me e utt (x, t) = c2 u(x, t) sur Rd R + u(x, 0) = f (x) pour x Rd ut (x, 0) = g(x) pour x Rd On dnit la moyenne sphrique de u : Mu (x, r, t) = e e avec u(x, t) = Mu (x, 0, t) . Grce a lquation de Darboux (2.15), on a a ` e or, par un calcul direct on obtient x Mu = 1 d
||=1
(2.16)
u(x + r, t) dS ,
2 d1 + 2 r r r
Mu
1 2 Mu x Mu = 2 . c t2
Donc Mu est solution, pour x Rd x, de lquation dEuler-Darboux-Poisson : e e 2 2 Mu (x, r, t) = c2 + d 1 M (x, r, t) pour (r, t) R R u + t2 r 2 r r (2.17) pour r R Mu (x, r, 0) = Mf (x, r) Mu (x, r, 0) = Mg (x, r) pour r R t On a donc transform lquation des ondes dans Rd en une quation aux drives partielles e e e e e avec variables scalaires r et t
39
2.4.4
On se propose de rsoudre le probl`me (2.16) pour d = 3 . e e En utilisant (2.17) on montre que la fonction v(r, t) = rMu (x, r, t) est solution de lquae tion des ondes pour d = 1 vtt (r, t) = c2 vrr (r, t) pour (r, t) R R + . v(r, 0) = rMf (x, r) pour r R vt (r, 0) = rMg (x, r) pour r R 1 v(r, t) r 1 1 = (ct + r)Mf (x, ct + r) (ct r)Mf (x, ct r) + 2r 2cr
ct+r
Mg (x, ) d .
ctr
Finalement, en faisant tendre r vers 0, on obtient la reprsentation suivante dune solution e u de (2.16) pour d = 3. u(x, t) = t Mf (x, ct) + t Mg (x, ct) , t (2.18)
et la formule de Kirchhoff pour le probl`me (2.16) en dimension 3 e u(x, t) = 1 4c2 t g(y) dSy +
|yx|=ct
1 t 4c2 t
f (y) dSy
|yx|=ct
(2.19)
Toute solution du probl`me (2.16) en dimension 3 qui est de classe C2 sur R3 R+ , est e donne par la formule (2.19) et est donc unique. e Rciproquement, pour f C3 (R3 ) et g C2 (R3 ), la fonction dnie par la formule de e e 2 3 Kirchhoff (2.19), est de classe C sur R R+ et est solution de (2.16). Grce au changement de variables y = x + ct , lquation (2.19) devient a e u(x, t) = 1 4 t g(x + ct ) + f (x + ct ) + ct
||=1
f (x + ct ) . dS .
(2.20)
Remarques : 1) Les formules (2.18) et (2.20) montrent que lon perd de la rgularit. e e 1 2 1 En eet si f C (R) et g C (R), alors Mf C (R ) et Mg C (R2 ) et u est de classe C1 pour t > 0 . Donc u est moins rguli`re que sa valeur initiale u(x, 0) ou, autrement dit, des petites e e irrgularits dans la condition initiale f peuvent saccumuler en t > 0 et rendre u(., t) e e moins rguli`re. e e Il est intressant dtudier le comportement de u par rapport a une nergie e e ` e E(t) = 1 2
R3
40 On montre que
1 dE , on a (t) = 0 pour tout t > 0 . 3 1 + |x| dt En particulier, si u = 0 pour |x| grand, le comportement L2 ne change pas avec t. Si on suppose que |ut (x, t) u(x, t)| 2) Grce a (2.20) on constate que la valeur de u(x, t) dpend des valeurs de f et g pour a ` e 2 y S (x, ct). Donc pour d = 3 le domaine de dpendance est la surface S 2 (x, ct). e
(x, t) t (y, t) ct R3 {t}
dE = dt
R3
ut (utt u) + c2 div(ut u) dx .
PSfrag replacements
Rciproquement, les donnes initiales au point (y, 0) R3 R+ nont de linuence sur e e u a linstant t quaux points x R3 qui vrient |y x| = ct. ` e Le cne dinuence de (y, 0) est la surface dans lespace-temps R3 R+ , forme par o e le cne de sommet (y, 0) et daxe de symtrie {(y, t)/t 0}, cest-`-dire lensemble o e a (x, t) / |y x| = ct, t > 0 . 3) Supposons que le support de f et g est = B 3 (0, r) R3 . Pour que u(x, t) soit non nul il faut que x S 2 (y, ct) avec y . Or lunion de toutes les sph`res de centre y et de rayon ct est la sph`re de centre e e 0, de rayon ct et dpaisseur 2r. e Donc supp(u) = B 3 (0, ct + r) \ B3 (0, ct r) et on constate que le support de u , dans R3 , se rpand a la vitesse nie c. e ` Pour x x et pour tout t < t1 = (|x| r)/c et t > t2 = (|x| + r)/c on a u(x, t) = 0 : le e front donde arrive en x a linstant t1 et a partir de linstant t2 londe est passe. ` ` e Dans R3 on a le Principe (fort) de Huygens : le domaine de dpendance est une surface dans R3 . e Cette proprit permet la modlisation par lquation des ondes de la transmission de ee e e 3 signaux a vitesse nie dans R . `
41
ct r
r
x1 x x
PSfrag replacements ct + r
4) Pendant que le support de u se rpand u dcroit en 1/t. e e En eet u(x, t) dpend des valeurs de f et g sur S 2 (x, ct) B 3 (0, r). Comme laire e maximale de cette intersection est laire de S 2 (0, r) = 4r 2 , et en supposant que g, f et ses drives partielles dordre 1 sont borns, on a grce a (2.20) e e e a ` |u(x, t)| r2 1 t sup |g| + sup |f | + ct sup | f | pour t grand. 2 t2 c t
2.4.5
Pour obtenir la solution de lquation des ondes pour d = 2 nous allons utiliser la mthode e e de descente de Hadamard : La solution de lquation aux drives partielles dans R2 est considre comme tant une e e e ee e 3 solution de lquation des ondes dans R indpendante de la variable x3 . e e La solution u(x1 , x2 , t) est donne par la formule de Kirchhoff (2.19) avec x3 = 0 . e Les conditions initiales sont f (y) = f (y1 , y2 ) et g(y) = g(y1 , y2 ) . Pour x = (x1 , x2 , 0) , y = (y1 , y2 , y3 ) on a ct = | x| = y Lintgrale sur les points y de la sph`re S 2 (, ct) devient une intgrale sur les y = (y1 , y2 ) e e x e du disque B2 (x, ct) , o` x = (x1 , x2 ). u On obtient la formule de Poisson pour lquation des ondes en dimension deux e 1 2c g(y) 1 dy + 2 t 2c |y x| f (y) dy . |y x|2 (y1 x1 )2 + (y2 x2 )2 + y3 2 .
u(x, t) =
|yx|<ct
c 2 t2
|yx|<ct
c 2 t2
42
Remarques : 1) Le domaine de dpendance de (x, t) est le disque B2 (x, ct). e Rciproquement, les donnes initiales en un point (y, 0) R2 R+ ont de linuence a e e ` linstant t sur tous les points de B2 (y, ct) . Le cne dinuence est le cne solide dans R2 R+ : (x, t) / |y x| ct, t > 0 . o o
t (y, t) ct R2 {t}
PSfrag replacements
(x, t)
Les perturbations en un point x R2 vont durer indniment : dans R2 on ne peut pas e transmettre de linformation. 2) Comparer le comportement des solutions de lquation de la chaleur et de lquation e e des ondes (d = 1, 2 et 3).
2.4.6
On sintresse au probl`me e e
Grce a la linarit il sut dtudier le probl`me avec conditions initiales homog`nes a ` e e e e e utt (x, t) c2 u(x, t) = w(x, t) sur Rd R + u(x, 0) = 0 pour x Rd ut (x, 0) = 0 pour x Rd
(2.21)
On a la
u(x, t) =
0
U (x, t, s) ds
pour x Rd et t 0,
43
Or on peut calculer V en dimension un, deux et trois : 1 solution de (2.21) pour d = 1 ; on a V (x, t, s) = 2c 1 2c
t
La mthode de Duhamel permet de rduire le probl`me (2.21) en une succession de proe e e bl`mes homog`nes. e e Vtt = c2 V pour x Rd , t > 0 Si on pose V (x, t, s) = U (x, t + s, s) on a V (x, 0, s) = 0 pour x Rd , t = 0 . Vt (x, 0, s) = w(x, s) pour x Rd , t = 0
x=ct
o` U (x, t, s) est, pour tout s 0, solution de u Utt = c2 U pour x Rd , t > s 0 U (x, s, s) = 0 pour x Rd , t = s Ut (x, s, s) = w(x, s) pour x Rd , t = s .
w(y, s) dy,
xct
do` u
u(x, t) =
w(y, s) dyds ;
0 |yx|<c(ts)
1 2c
|yx|<ct
c 2 t2
w(y, s) , dy , |y x|2
do` u
u(x, t) =
on obtient
|yx|<ct
w y, t 1 |y x| c dy (potentiel retard). e |y x|
44
2.5
Pour u C2 (), ouvert de Rd , une quation aux drives partielles linaire dordre 2 e e e e scrit en tout point x de : e 2u aij (x) (x) + xi xj 1i,jd
d
bi (x)
i=1
(2.22)
o` lon peut supposer aij = aji , pour tout 1 i, j d. u On dnit le symbole de lquation aux drives partielles (2.22) , en x et pour tout e e e e d = (1 , . . . , d ) R , par aij (x) i j . x () =
1i,jd
Pour tout x , x est une forme quadratique de matrice associe A(x) = aij (x) e Considrons un diomorphisme : avec x = (x) et dt(D) = e e e Pour x x, on pose P = Dx t , i.e. pij = e 2u (x ) = xi xj =
1k,ld
1i,jd
xi xj
= 0.
1i,jd
1i,jd
aij (x )
1k,ld
1i,jd
akl ( ) x
2u ( ) + (termes dordre 1) . x xk xl
Cest-`-dire que les coecients dordre deux de lquation aux drives partielles (2.22) a e e e sont transforms de la mme faon que les coecients de la forme quadratique x sous e e c la transformation = P : on a x () = t P t A(x )P . Or, grce a un changement de variables adapt, une forme quadratique relle peut tre a ` e e e rduite a sa forme canonique : e `
m r 2 i i=1
2 i , i=m+1
o` 1 m r d. Donc, en un point x donn, lquation aux drives partielles u e e e e (2.22) peut toujours scrire sous la forme e
m
i=1
2u 2u ( ) x ( ) + x x2 i x2 i i=m+1
x i ( ) u ( ) + c( ) u( ) = d( ) . x x x b x xi i=1
45
Dfinition 2.5.1 e On dit que lquation aux drives partielles (2.22) est : e e e e elliptique au point x si r = d et m = 0 ou m = d, i.e. A(x ) est dnie positive ou ngative et a d valeurs propres relles non nulles de mme signe ; e e e hyperbolique au point x si r = d et m = 1 ou m = d 1, i.e. A(x ) est indnie et a e d valeurs propres relles non nulles dont une de signe contraire aux d 1 autres ; e parabolique au point x si r = d 1 et m = 0 ou m = d 1, i.e. A(x ) a 0 comme valeur propre et d 1 valeurs propres relles non nulles de mme signe. e e Exemples : 1. Lquation de Laplace dans Rd est elliptique en tout point : e
d
x Rd , Rd : x () =
i2 .
i=1
(x, t) Rd R , Rd+1 : (x,t) () = + 3. Lquation des ondes est hyperbolique en tout point : e
d
i2 .
i=1
i=1
2 i2 d+1 .
Elle est donc elliptique si x2 > 0, hyperbolique si x2 < 0 et parabolique si x2 = 0 . 5. Lquation ux1 x2 = 0 est hyperbolique dans R2 . e Remarques : 1. Pour tout point x Rd on a une transformation qui rduit (2.22) sous forme canoe nique or cette transformation dpend de x . Pour d 3 il est en gnral impossible e e e de trouver un diomorphisme qui permet dcrire lquation aux drives pare e e e e tielles sous forme rduite sur une rgion donne (arbitrairement petite). e e e Dans le cas d = 2 et sous quelques hypoth`ses sur les coecients de (2.22) on peut e trouver qui permet de mettre (2.22) sous forme rduite. e 2. La classication ci-dessus a t donne pour des quations aux drives partielles ee e e e e linaires dordre deux. On retrouve ces mmes dnominations pour des quations e e e e aux drives partielles plus gnrales, mais qui partagent les proprits essentielles e e e e ee avec les quations aux drives partielles linaires prsentes dans ce chapitre. e e e e e e
46
On largit la classe des quations aux drives partielles hyperboliques aux quae e e e e tions qui admettent des solutions de type ondes, i.e. on a une vitesse de propagation nie. Considrons par exemple le syst`me dquations aux drives partielles e e e e e dordre un en t Ut (x, t) + A Ux (x, t) = F (x, t) dans R R , + o` U = ( u1 uN )t est dni sur R R+ et A est une matrice constante u e dordre N N . On a un syst`me hyberbolique a coecients constants si A est e ` diagonisable sur R . Cette dnition se gnralise au cas o` A(x, t) et pour U : Rd R+ RN . e e e u Lquation de transport est une quation aux drives partielles hyperbolique e e e e dordre un en t . Pour les quations aux drives partielles elliptiques a coecients non constants e e e ` on utilisera un crit`re uniforme, voir page 61 . e
Dans ce chapitre on sintresse au probl`me suivant : e e d Soit un ouvert de R , on veut trouver u : Rd R solution de F x, u(x), u(x) = 0 u(x) = h(x) pour x pour x . (3.1)
On suppose que F et h sont au moins de classe C2 et que est une hypersurface rguli`re. e e c Rd+1 , on obtient lquae 1 tion de transport, cf. page 15, que lon int`gre sur les droites caractristiques. e e Exemple : Pour F (x, u, u) = v. u f = 0, avec v = Ide : soit u une solution de (3.1) et x , comme pour lquation de transport, on e e veut calculer u(x) a partir de u(x0 ) = h(x0 ) , x0 , en intgrant le long dun chemin ` e X(s) .
PSfrag replacements u h(x0 ) u(x) z(s)
X(s) X(0) = x0
48 et
u(X(s)) .
Sur la courbe X lquation (3.1) devient : F x1 (s), . . . , xd (s), z(s), g1 (s), . . . , gd (s) = 0 . e Fg 1 F x1 . . On note X F = . et G F = . . . . F xd Fg d
3.2
Pour une solution u de (3.1), le thor`me qui suit fait le lien entre la courbe paramtre e e e e X, la valeur de u sur X, z, et la valeur du gradient de u le long de X, G. Thor`me 3.2.1 (Syst`me caractristique) e e e e 2 Soit u C () une solution de (3.1). Si X est solution de X (s) =
G F (X(s), z(s), G(s))
(3.2.1)
alors, pour tout s tel que X(s) , z et G sont solution de z (s) = G (s) = Remarques : 1. Les 2d+1 quations direntielles ordinaires (3.2.1-3) forment le syst`me dquations e e e e caractristiques associ a lquation aux drives partielles dordre un (3.1). e e` e e e 2. Les fonctions X, z et G sont appeles les caractristiques et X est la projection des e e caractristiques ou la courbe caractristique. e e Si lon se donne des conditions initiales (X(0), z(0), G(0)) , alors, grce a la rgularit de a ` e e F , on a existence dune solution locale unique (X, z, G) du syst`me autonome dquations e e direntielles ordinaires caractristiques. e e Donc, pour obtenir u, solution de (3.1), il faut rsoudre (3.2.1-3) pour X(0) = x0 . e En un premier temps, pour simplier ltude, on va montrer que lon peut supposer que e est plat au voisinage de x0 . En eet, comme est rgulier il est possible daplatir au voisinage de x0 , cf. annexe e B. Il existe , diomorphisme de Rd , et > 0 tel que pour tout y Bd (x0 , ) : e (y) = (y1 , . . . , yd1 , 0). Si on pose x = (x) et u() = u(1 ()), alors u vrie lquation aux drives partielles x x e e e e dordre un x F (, u, u) = 0 pour x () u = h pour x () . Donc, si est plat au voisinage de x0 , les conditions initiales en x0 sont : X(0) = x0 , z(0) = h(x0 ) = h(x0 , . . . , x0 , 0) et G(0) = p0 . 1 d1
G F (X(s), z(s), G(s)).G(s) X F (X(s), z(s), G(s))
(3.2.2) (3.2.3)
49
Or pour x0 donn, z(0) = h(x0 ) est unique, il reste a trouver p0 . e ` En drivant u(y1 , . . . , yd1 , 0) = h(y1 , . . . , yd1 ) en x0 et en notant p0 = (p0 , . . . , p0 ) = e 1 d G(0) , on obtient p0 = hxi (x0 ) , pour i = 1, . . . , d 1 i F (x0 , h(x0 ), p0 , . . . , p0 ) = 0 1 d On dit que (x0 , h(x0 ), p0 ) est admissible si p0 = (p0 , . . . , p0 ) est solution de (3.3). 1 d Or on veut rsoudre (3.1) dans un voisinage de x0 . e Donc, pour tout y = (y1 , . . . , yd1 , 0) Bd (x0 , ), on va rsoudre (3.2.1-3) avec les e conditions initiales X(0) = y , z(0) = h(y) et G(0) = q(y). Il faut donc trouver q : Bd (x, ) Rd telle que q(x0 ) = p0 , o` p0 vrie (3.3). u e De plus il faut assurer que (y, h(y), q(y)) est admissible, do` u qi (y) = hxi (y) , pour i = 1, . . . , d 1 F (y, h(y), q1(y), . . . , qd (y)) = 0 et on a le lemme suivant qui assure lexistence de q : Lemme 3.2.1 Soit (x0 , h(x0 ), p0 ) admissible et supposons que Fgd (x0 , h(x0 ), p0 ) = 0 . Alors il existe une solution unique q de (3.4) et (3.5). e e La condition Fgd (x0 , h(x0 ), p0 ) = Xd (0) = 0 , du lemme prcdent signie que le champ de vecteurs X est transverse a en x0 . ` Si nest pas plat au voisinage de x0 cette condition devient
G F (x0 , h(x0 ), p0 ).n(x0 )
(3.3)
(3.4)
(3.5)
= 0,
o` n(x0 ) est la normale a en x0 . On dit que est non caractristique en x0 pour (3.1). u ` e Dans ce cas X (0).n(x0 ) = 0 , cest-`-dire X (0) Tx0 . a
Rd PSfrag replacements
y n(x0 ) x0
X (0)
Finalement, on montre quen recollant les solutions de (3.2.1-3) pour des conditions initiales admissibles, on obtient une solution locale unique de (3.1) et on a le
50
Thor`me 3.2.2 (Existence locale) e e Soit x0 , (x0 , h(x0 ), p0 ) admissible et G F (x0 , h(x0 ), p0 ).n(x0 ) = 0 . Alors il existe un voisinage V de x0 dans Rd et une fonction unique u C2 (V ), tels que F (x, u(x), u(x)) = 0 u(x) = h(x) pour x V . pour x V .
3.3
quation aux drives partielles linaire dordre un : e e e e On a la forme gnrale e e F (x, u, u) = a(x). u(x) + b(x)u(x) + c(x) = 0 .
Comme G F = a(x) , est non caractristique en x0 si a(x0 ).n(x0 ) = 0 . e Les quations direntielles ordinaires (3.2.1) et (3.2.2) donnent le syst`me caractristique e e e e rduit e X (s) = a(X(s)) z (s) = b(X(s)) z(s) c(X(s)) . (3.6)
Ce syst`me direntiel est indpendant de G et grce au rsultat dunicit dune solution e e e a e e locale rguli`re il sut dintgrer (3.6) pour obtenir cette solution. e e e Exemple : On consid`re lquation aux drives partielles e e e e x1 u x 2 x 2 u x 1 = u u(x) = h(x) dans (R )2 + . pour x = R {0} . +
On montre que pour tout x0 > 0, a(x0 , 0).n(x0 , 0) = x0 > 0, donc est partout non caractristique. Les courbes caractristiques sont des quarts de cercles et la solution du e e probl`me est e (x1 , x2 ) (R )2 + : u(x1 , x2 ) = h x2 1 + x2 2 e
arctan
x2 x1
quation aux drives partielles quasilinaire dordre un : e e e e La forme gnrale est e e F (x, u, u) = a(x, u). u(x) + b(x, u) = 0 .
On a G F = a(x, u) et est non caractristique en x0 si a(x0 , h(x0 )).n(x0 ) = 0 . e Le syst`me caractristique rduit scrit e e e e X (s) = a(X(s), z(s)) z (s) = b(X(s), z(s)) .
3.3. CAS LINEAIRE ET QUASILINEAIRE
(3.7)
51
Comme dans le cas linaire, il sut de rsoudre (3.7) pour trouver la solution rguli`re e e e e locale. Exemple : On consid`re lquation aux drives partielles e e e e ux 1 + u x 2 = u 2 u=h dans = R R + . sur = R {0} .
On montre que est non caractristique en tout point et que les courbes caractristiques e e sont des droites. La solution est donne par e u(x1 , x2 ) = h(x1 x2 ) 1 x2 h(x1 x2 ) pour tout (x1 , x2 ) tant que x2 h(x1 x2 ) = 1 .
Dans la suite on va tudier dautres cas o` la solution nexiste pas dans tout . e u
3.4
Dans cette section on sintresse a un type particulier dquations aux drives partielles e ` e e e quasilinaires dordre un, les lois de conservation scalaires. Ce sont des probl`mes de la e e forme ut (x, t) + f (u(x, t))x = 0 o` f C2 (R) et avec u(x, 0) = h(x) . u sur R R + (3.8)
3.4.1
Modlisation e
On consid`re des particules en dplacement le long de laxe Ox ; u(x, t) est la densit de e e e particules en x a linstant t. Soit =]a, b[, on pose `
b
m (t) =
a
u(x, t)dx .
On suppose que la variation de masse dans est uniquement due au ux de particules, f (u), en x = a et x = b . En tenant compte du signe du ux on obtient : d m (t) = dt
b
ut (x, t)dx
a
f (u).n d
f (u)x dx
a
div(f (u)) dx .
En supposant que u est rgulier et comme est quelconque, on obtient la loi de consere vation de la masse ut + f (u)x = 0 . Si lon se donne la vitesse v des particules alors f (u) = vu ([p/s] = [p/m][m/s]).
52 Exemples :
1. Si la vitesse est constante, v = c, on retrouve lquation de transport. e 2. Si on modlise le trac routier, une relation simple entre la densit du trac u et la e e vitesse des voitures (ou vlos) v est e v = vmax 1 u umax ,
Pour une route vide (u = 0) la vitesse est vmax mais la vitesse dcroit vers 0 quand e le trac arrive a saturation (u = umax ), le ux est ` f (u) = u vmax 1 u umax .
Des expriences amricaines, menes en 1959 dans un tunnel, ont conclu que le ux e e e de vhicules y pouvait tre modlis par f (u) = au ln(umax /u) . e e e e
3.4.2
Modlisation : e On consid`re des particules en dplacement le long de laxe Ox, v(x, t) est la vitesse de la e e particule qui se trouve en x a linstant t. On se donne la vitesse initiale h(x) = v(x, 0). ` On consid`re une particule X qui se trouve en x0 = X(0) pour t = 0 et en x = X(t) e a linstant t, alors v(X(t), t) est la vitesse de la particule X a linstant t. En particulier ` ` v(x0 , 0) = v(X(0), 0) = h(x0 ) . Supposons que chaque particule se dplace a vitesse constante : v(X(t), t) = h(X(0)) . e ` d En drivant par rapport a t (` particule xe ), on a : vx X(t) + vt = vx v + vt = 0 . e ` a dt On obtient ainsi lquation de Burgers pour des uides non visqueux : e vt + vvx = 0 dans R R + v(x, 0) = h(x) pour x R Cest une quation du type (3.8) avec f (v) = 1 v 2 . e 2 Etude de lquation : e = R {t = 0} et on montre que est non caractristique en tout point (x, 0). Pour tout e x R il existe un voisinage de (x, 0) R2 dans lequel (3.9) admet une solution unique rguli`re. e e x (s) = z x(0) = x0 t (s) = 1 avec les c.i. t(0) = 0 Le syst`me caractristique scrit : e e e . z (s) = 0 z(0) = h(x0 ) (3.9)
Les courbes caractristiques sont les droites passant par (x0 , 0) dquation e e x = x0 + h(x0 ) t et le long desquelles v(x, t) = h(x0 ). Une droite caractristique reprsente le chemin dans lespace-temps de la particule X, e e partie de x0 = X(0) .
53
Comme v(x, t) = h(x0 ) = h(x h(x0 ) t) = v(x v(x, t) t), la solution de (3.9) vrie e lquation implicite e v = h(x v t) . Supposons quune deuxi`me particule X1 , X1 (0) = x1 > x0 , se trouve a linstant t en x : e ` x = x0 + h(x0 ) t = x1 + h(x1 ) t donc si h(x0 ) = h(x1 ) : t= x1 x 0 . h(x1 ) h(x0 )
Si h est croissante il nexiste pas de t > 0 vriant cette relation. e Si h nest pas croissante un tel instant t > 0 existe et v nest plus une fonction, v(x, t) devant tre gal a h(x0 ) et h(x1 ). e e `
t t PSfrag replacements
x0
x1 x
x0
x1 h croissant
h non croissant
Les particules rapides (vitesse h(x0 )) vont rattraper les particules lentes (vitesse h(x1 )).
PSfrag replacements
v(x, 0)
v(x, t)
X(0) = x0 X1 (0) = x1
x = X(t) = X1 (t)
Calculons vx le long dune droite caractristique : vx = e Donc si h (x0 ) < 0 pour t <
h (x0 ) . 1 + h (x0 )t
1 on a |vx | qui tend vers linni. h (x0 ) 1 1 = Le premier instant o` v nest plus dni correspond a T = min u e ` , xR h (x) h (x ) o` x est un minimum de h . On dit que la solution explose (angl. blow-up) en T . u Il ne peut y avoir de solution C1 au del` de linstant T . a Pour pouvoir dnir des solutions au del` de T , i.e. sur R R , il faut introduire des e a + solutions gnralises. e e e
54
3.4.3
Soit la fonction test D(R R+ ), en supposant u de classe C1 et en multipliant (3.8) par on obtient par intgration par parties e
+
0=
0 R
ut + f (u)x dxdt
+ +
=
R
u(x, 0)(x, 0) dx
ut dtdx
f (u)x dxdt
0 R
Cette derni`re relation est valable si u nest pas C1 , do` la e u Dfinition 3.4.1 e On dit que u L (R R+ ) est une solution faible du probl`me e ut + f (u)x = 0 dans R R + u(x, 0) = h(x) pour x R si pour toute fonction test D(R R+ ) on a
+ 0 R
(3.10)
ut + f (u)x dxdt =
h(x)(x, 0) dx
R
(3.11)
Nous allons tudier les solutions gnralises pour un cas particulier simple. e e e e Soit u une solution gnralise de (3.10) dans un ouvert R R , et C une courbe e e e + rguli`re qui partage en a gauche et + a droite de C . e e ` ` On va supposer que u est de classe C1 sur et sur + , C est lensemble de discontinuit e de u .
t C
PSfrag replacements +
En choisissant les fonctions test a support compact dans , resp. + , on montre que u ` est une solution forte de (3.10) dans , resp. + . Pour une fonction test dont le support est partag en deux par C on trouve en utilisant e la relation (3.11) (f (u ) f (u+ ))n1 + (u u+ )n2 d = 0 ,
n1 est le vecteur normal unit extrieur a et u , resp. u+ , est la limite de e e ` n2 u a gauche, resp. a droite, de C . ` ` o` n = u
55
qui fait le lien entre la vitesse de propagation de la discontinuit et la valeur du saut de e u, resp. f (u), en ((t), t)). On voit sur des exemples que la condition de Rankine-Hugoniot permet en gnral e e plusieurs solutions gnralises. Nous allons introduire une contrainte supplmentaire pour e e e e liminer des solutions non physiques . e x (s) = f (z) t (s) = 1 Le syst`me caractristique de (3.10) scrit e e e . z (s) = 0 On en dduit que u(x, t) = h(x0 ) le long des droites caractristiques x = x0 + f (h(x0 )) t, e e tant que u est de classe C1 . Considrons c C avec u (c) = u+ (c) et tel que c est le premier point de rencontre dune e droite caractristique venant de (x0 , 0) avec une droite caractristique venant de e e (x1 , 0) + . Alors u (c) = h(x0 ) et u+ (c) = h(x1 ).
t PSfrag replacements c + C
x0
x1
i.e. a chaque instant la vitesse du ux a gauche de la discontinuit est plus grande qu` ` ` e a droite. Dfinition 3.4.2 e Une discontinuit C est un choc si les conditions (RH) et (E) sont vries. e e e
56
3.4.4
On reprend lexemple du trac routier o` u(x, t) est la densit de voitures et vlos, et u e e v(x, t) la vitesse en x a linstant t. La distribution initiale de vhicules est donne par ` e e u(x, 0) = h(x) , o` h : R [0, 1] et il ny a ni bretelles dacc`s, ni bretelles de sortie sur u e cette route. En supposant que v = 1 u, i.e. umax = vmax = 1, on a f (u) = u(1 u), f (u) = 1 2u et on obtient ut + (u(1 u))x = 0 u(x, 0) = h(x) pour (x, t) R R + pour x R . (3.14)
Les droites caractristiques ont pour quation : x = x0 + (1 2h(x0 ))t . e e Si h est de classe C1 on sait quil existe une solution unique u C1 (R]0, T [), pour 1 . T = inf x0 R 2h (x0 ) Si h est dcroissante T = + , i.e. initialement il y a plus de vhicules a gauche qu` e e ` a droite et les vhicules a gauche vont moins vite que ceux a droite et ne pourront jamais e ` ` les rattraper. Si h nest pas dcroissante a linstant T la solution C1 cesse dexister : des vhicules a e ` e ` grande vitesse arrivent en une zone a grande densit et faible vitesse. ` e Considrons le probl`me (3.14) avec la condition initiale h(x0 ) = u(x0 , 0) = e e un probl`me de ce type est appel probl`me de Riemann. e e e En utilisant f et f on montre que la condition (RH) scrit = 1 (u + u+ ) , e la condition (E) devient ici u < u+ . Si h(x0 ) = t 1/2 x0 < 0 , alors (t) = et u(x, t) = 1 x0 > 0 2 1/2 x < t/2 . 1 x > t/2 h G x0 < 0 h D x0 > 0
` A droite de zro il y a un bouchon (v = 0) a linstant t = 0 et les vhicules venant de la e ` e gauche simmobilisent au fur et a mesure que londe de choc se dplace vers la gauche. ` e Sur la gure suivante on a reprsent, a gauche, les droites caractristiques. e e ` e
t t
PSfrag replacements
0 Caractristiques e
x0
Sur la gure de droite sont reprsents les trajectoires des vhicules, obtenues comme e e e
57
suit : Notons X(t) la position a linstant t dun vhicule qui est parti en X(0) = x0 < 0 . ` e dX t/2 + x0 0 t < x0 On conna t = v(x, t) , en intgrant on obtient X(t) = e x0 /2 t > x0 dt reprsent dans le plan xt . e e Si h(x0 ) = 1/6 x0 < 0 , on trouve u(x, t) = 1/3 x0 > 0 1/6 x < t/2 . 1/3 x > t/2
Dans ce cas londe de choc se dplace vers la droite, les vhicules rapides (v = 5/6) venant e e de gauche doivent freiner pour rouler dans la zone dense a v = 2/3 . On reprsente sur la ` e gure suivante les droites caractristiques et les trajectoires des vhicules. e e
t t
PSfrag replacements
0 Caractristiques e
x0
Soit h(x0 ) =
1 x0 < 0 , si on cherche une solution de (3.14) avec une disconti1/2 x0 > 0 1 x < t/2 nuit vriant (RH) on trouve u(x, t) = e e . 1/2 x > t/2 Cette discontinuit nest pas un choc car (E) nest pas vri. Les vhicules a gauche de e e e e ` zro sont a larrt et dmarrent au fur et a mesure, leur vitesse passant instantanment e ` e e ` e de 0 a 1/2. On remarquera que les caractristiques sortent de la discontinuit . ` e e
t
PSfrag replacements
0 Caractristiques e
x0
58
Dans la zone {(x, t) / t > 0, t < x < 0} on a un phnom`ne de rarfaction ; pour pouvoir e e e acclrer les vhicules doivent avoir de lespace : la vitesse augmente de faon rguli`re de ee e c e e v = 0 a v = 1/2 et la densit diminue. ` e
1 x < t 1/2 x/2t t < x < 0 . Un solution physiquement plus acceptable est u(x, t) = 1/2 x>0
t t t0 x = t x = t
PSfrag replacements
0 Caractristiques e x x0
x0
f v dx
(*)
Sous les hypoth`ses de rgularit, cites ci dessus, () est quivalent a lquation aux e e e e e ` e drives partielles u = f . e e Or pour pouvoir crire () il sut que u, v, f , uxi et vxi soient dans L2 (). On obtient e ainsi une formulation faible de u = f . Il faudra cependant complter cette formulation e faible pour intgrer la condition de Dirichlet homog`ne. e e
4.1
Espaces de Sobolev
Dans toute cette section sera un ouvert de Rd . Dfinition 4.1.1 e On appelle espace de Sobolev sur , lespace fonctionnel W k,p() = u Lp () Zd , || k : D u Lp () ,
60
1 2
u. v dx
la norme associe. e
Rappelons que D() est dense dans L2 (). On montre quen gnral D() nest pas dense dans H 1 () e e et on introduit ladhrence de D() dans H 1 (), on note e
1 H0 () = D() H 1 ()
Exemples :
1 1. Pour = Rd on montre que H0 (Rd ) = H 1 (Rd ) .
2. Soit un ouvert born de Rd avec rgulier. Si u C1 () et u| = 0 alors u e e 1 H0 (), cest-`-dire quil existe une suite (n ) dans D() tel que lim u n H 1 () = 0. a
n 1 On interpr`te H0 () comme lespace des fonctions de H 1 () qui sannulent sur . On e verra plus loin quel sens prcis donner a cette armation. e `
Proposition 4.1.1 (Ingalit de Poincar) e e e Soit un ouvert born de Rd , alors il existe une constante C() telle que e
1 u H0 () :
u2 dx C()
| u|2 dx .
Remarque : Soit un ouvert born de Rd et u = I , pour R . e 1 1 Alors u H () pour tout mais par lingalit de Poincar on montre que u H0 () e e e si et seulement si = 0 , i.e. u est la fonction constante nulle sur . 1 On a donc en gnral H0 () H 1 () . e e = Proposition 4.1.2 Soit un ouvert born de Rd , alors e u
H 1 ()
(u + | u| )dx
1 2
et
1 H0 ()
| u| )dx
1 2
61
Thor`me 4.1.2 (Thor`me des traces) e e e e Soit un ouvert born de Rd avec rgulier. e e Alors il existe un oprateur linaire born T : H 1 () L2 () tel que e e e u C0 () H 1 () : T u = u| u H 1 () : T u L2 () C u o` C ne dpend que de . u e On dit que T u est la trace de u sur et on note T u = u| . Sous les hypoth`ses du thor`me des traces on a : e e e
1 u H0 () si et seulement si T u = 0 sur ; H 1 ()
pour tout u, v de H 1 () :
uxi v dx =
uvni dS
uvxi dx .
4.2
Dans la suite est un ouvert born de Rd et est rgulier. e e On consid`re le probl`me e e Lu = f sur u = 0 sur o` loprateur direntiel linaire L est dni par u e e e e
d
(4.1)
Lu =
xj
+
i=1
avec aij (x) = aji (x) et en supposant que les aij , bi , c sont dans L () et f L2 (). On suppose de plus que L est elliptique au sens suivant : > 0 , Rd : aij (x) i j ||2 p.p. .
1i,jd
Soit u une solution forte de (4.1) et v D(), alors par intgration par parties on obtient e
d
(Lu)v dx =
f v dx ,
f v dx .
1 Cette derni`re relation a un sens mme si lon remplace u et v par des fonctions de H0 () , e e do` la dnition : u e
62
a(u, v) =
f v dx ,
o` a(., .) est la forme bilinaire associ a loprateur elliptique L et dnie pour tout u, v u e e` e e 1 de H0 () par
d
a(u, v) =
1i,jd
(4.2)
On va prsenter maintenant des rsultats abstraits danalyse fonctionnelle qui vont sape e pliquer notre probl`me. e e Thor`me 4.2.1 (Theor`me de Lax-Milgram) e e Dans un espace de Hilbert rel H , dont la norme est not . , on consid`re une forme e e e bilinaire a(., .) vriant e e (i) > 0 , u, v H : (ii) > 0 , u H : |a(u, v)| u v u 2 a(u, u) .
Alors pour toute forme linaire continue, l H , il existe un unique u H vriant e e v H : a(u, v) = l(v) . (V)
On dit que (V ) est un probl`me variationnel. e Thor`me 4.2.2 e e Sous les hypoth`ses du thor`me de Lax-Milgram et en supposant que la forme bilinaire e e e e a(., .) est symtrique on montre que la solution u H de (V ) est lunique minimum de la e fonctionnelle 1 J[v] = a(v, v) l(v) , 2 cest-`-dire J[u] = min J[v] . a
vH
Application :
1 Le produit scalaire sur H0 (), a(u, v) =
1 1 Il existe donc un unique u H0 () vriant : v H0 () : a(u, v) = l(v) . e On a montr lexistence et lunicit dune solution faible de lquation de Poisson avec e e e conditions de Dirichlet homog`ne, e
u = f sur u = 0 sur
63
o` f L2 () est donne. u e De plus cette solution est le minimum de J[v] = cf. principe de Dirichlet.
1 | v|2 f v dx 2
Remarque : Si a(., .) est symtrique on na pas besoin du thor`me de Lax-Milgram, le e e e thor`me de reprsentation de Riesz sapplique directement. e e e Proposition 4.2.1 Soit a(., .) la forme bilinaire associ a loprateur elliptique L et dnie par (4.2), on a : e e` e e
1 1) > 0 , u, v H0 () :
2) > 0 , 0 , u
1 H0 ()
|a(u, v)| u
2 1 H0 ()
1 H0 ()
a(u, u) + u
1 H0 () 2 L2 ()
Cette proposition montre que pour L on a immdiatement la condition (i) du thor`me e e e de Lax-Milgram, par contre la condition (ii) nest pas vrie en gnral. e e e e On a le thor`me : e e Thor`me 4.2.3 e e Il existe 0 tel que pour tout et pour tout f L2 (), il existe une solution 1 faible unique u H0 () du probl`me e Lu + u = f sur u = 0 sur .
4.3
Soit un ouvert born de R2 tel que est une union nie de polygones. e Soit Th un ensemble ni de polygones K, dintrieurs non vide et tels que e = K
KTh
Si K1 = K2 : int(K1 ) int(K2 ) = et K1 K2 est vide ou rduit a un sommet commun e ` ou cest un ct commun a K1 et K2 . oe ` On dit alors que Th est une triangulation de et h = max diam(K) .
KTh
On note Q1 = { P (x1 , x2 ) = a + bx1 + cx2 + dx1 x2 / a, b, c, d R } , lensemble des polynmes sur R2 de degr infrieur a un par rapport a chaque variable. o e e ` ` Cest un espace vectoriel de dimension 4, en particulier la valeur de P Q1 aux quatre sommets dun rectangle K dtermine enti`rement P sur K . e e
On note h lensemble des sommets des polygones de Th . Pour simplier on suppose que tous les polygones de la triangulation Th sont des rectangles.
64
Exemples : On prend K = [0, 1]2 . Ecrire lquation des polynmes suivants : e o 1. soit P1 Q1 tel que P (0, 0) = P (1, 1) = 1 et P (0, 1) = P (1, 0) = 0 ;
PSfrag replacements
x2
(0, 0) x1 (1, 0)
(0, 0) (1, 0)
x1
65
En eet, pour i {1, 2} la drive au sens des distributions de v par rapport a xi est la e e ` 2 fonction vi L () , dnie sur par vi = v|K xi pour tout K Th . e Soit h = {s1 , . . . , sN } , lensemble des sommets intrieurs a , on introduit les e ` fonctions i Vh telles que i (sj ) = ij . On montre que {i }1iN est une base de Vh .
N
Pour tout v Vh :
v=
i=1
v(si ) i .
(0, 0, 1) (s) =
5
1 si s = (0, 0) 0 sinon
0.
(1, 1)
(0, 1)
(1,1)
En utilisant le thor`me de Lax-Milgram on obtient la e e Proposition 4.3.2 Soit Vh un sous espace vectoriel de dimension nie dun espace de Hilbert rel H , e soit a(., .) une forme bilinaire dnie sur H vriant les conditions (i) et (ii) du thor`me e e e e e de Lax-Milgram et soit l H . Alors il existe u H et uh Vh solutions uniques des probl`mes e v H : vh Vh : et telles que u uh a(u, v) = l(v) a(uh , vh ) = l(vh ) inf u vh . vh Vh
1 Cette proposition sapplique immdiatement au cas o` H = H0 () et pour Vh dni dans e u e la proposition 4.3.1. Dans ce cas on obtient en plus la
Proposition 4.3.3 Soit D() , on dnit la projection, h , de sur Vh par : si h , h (si ) = (si ) . e Alors 1 lim h H0 () = 0 .
h0
66
u = f sur u = 0 sur
avec f L2 () donn. e 1 1 Comme D() est dense dans H0 () il existe D() tel que u H0 () < /2 pour > 0 quelconque. 1 De plus h Vh et, pour h susamment petit, h H0 () < /2. Donc u uh cest-`-dire a
1 H0 ()
u h
1 H0 ()
1 H0 ()
+ h
1 H0 ()
<,
h0
lim u uh
1 H0 ()
=0.
Calcul de uh :
N
uh Vh donc uh =
vh Vh : j = 1, . . . , N :
N
a(uh , vh ) =
f vh dx f j dx = Fj
a(uh , j ) = a(i , j ) Ui = Fj
i=1
j = 1, . . . , N :
Si on note A = (Aij )1i,jN = (a(j , i ))1i,jN , F = ( F1 . . . FN )t et U = ( U1 on est amen a rsoudre un syst`me linaire de taille N : e` e e e AU = F .
. . . UN )t
Note : plus h est petit, plus lapproximation de u est bonne mais le nombre de sommets augmente et donc aussi la taille du syst`me linaire. e e Algorithme 1) 2) 3) 4) Construction de la triangulation du domaine . Evaluation des i (1 i N ) . Calcul de Fj (1 j N ) par des mthodes de quadrature en dimension d = 2 . e Calcul de i j i j Aij = + dx1 dx2 . x1 x2 K x1 x2 KT
h
Or si si et sj ne sont pas sommets dun mme K alors supp(i ) supp(j ) = , la e matrice A est creuse. On peut ordonner les sommets de faon a obtenir une matrice bande. c ` 5) Dterminer U en rsolvant A U = F en utilisant des mthodes de rsolution nume e e e e riques (directes ou itratives). e 6) Pour augmenter la prcision raner Th (i.e. diminuer h) et recommencer en 1 . e
Dans ce chapitre nous allons tudier la mthode des dirences nies pour rsoudre nue e e e mriquement des quations aux drives partielles dvolution linaires dordre 1 en t . e e e e e e Le reprsentant le plus simple qui nous servira dexemple type pour illustrer les dnie e tions et techniques proposes, est lquation hyperbolique dordre un en t : lquation de e e e transport. Linconnue est une fonction u dnie pour (x, t) R R+ . Le domaine de dnition de e e u est discrtis par la grille Gh,k = {(xm , tn ) / xm = mh, m Z, tn = nk, n N} e e avec h et k des rels strictement positifs que lon choisira petits. e On dit que h est le pas de discrtisation spatiale et k le pas de temps . e La fonction u qui est dnie pour les variables continues (x, t) , prend la valeur e un = u(mh, nk) au point (xm , tn ) = (mh, nk) de la grille Gh,k . m Pour bien distinguer la solution continue u du rsultat du schma numrique, on va noter e e e n v la solution numrique qui est dnie uniquement sur Gh,k : (vm )mZ,nN . e e n Donc vm est la solution approche, au point (xm , tn ) , de lquation aux drives partielles. e e e e Pour obtenir un probl`me discret on remplace les drives partielles par des dirences e e e e nies, ainsi : u (xm , tn ) x u (xm , tn ) x u (xm , tn ) x 2u (xm , tn ) x2 un u n m+1 m , h un u n m m1 , h un u n m+1 m1 , 2h dirence nie progressive ; e dirence nie rtrograde ; e e dirence nie centre ; e e
Ces approximations sont obtenus grce a la formule de Taylor. On fait de mme pour a ` e les drives partielles par rapport a t : e e ` u (xm , tn ) t un+1 un m m , k dirence nie progressive en t . e
68
Exemple : Considrons lquation de transport ut + cux = 0 dans R R , avec la condition initiale e e + u(x, 0) = (x) pour x R et on suppose que c > 0. On va discrtiser cette quation aux drives partielles en utilisant un schma progressif e e e e e en espace et en temps :
n n+1 n v n vm vm v m + c m+1 =0. k h
(5.1)
n+1 n n n En posant = k/h on obtient : vm = (1 + c)vm cvm+1 = (1 + c cTh ) vm , o` T est loprateur de translation spatiale dni par T u(x, t) = u(x , t) pour tout u e e n n R , i.e. Th vm = vm+1 . On a par rcurrence : e 0 n vm = (1 + c c Th )n vm n
=
p=0 n
=
p=0
Ainsi le domaine de dpendance de v au point (xm , tn ) = (ih, nk) est form par les points e e xm , xm + h, xm + 2h, . . . , xm + nh . Or, on sait que u(x, t) = (x ct) et que le domaine de dpendance de un = u(xm , tn ) est restreint au point xm cnh . Le schma discret ne e e m tient donc pas compte des proprits de la solution de lquation aux drives partielles ee e e e n et la solution discr`te vm est en gnral dirente de un : le schma (5.1) ne peut pas e e e e e m tre convergent. e
t tn+1 PSfrag replacements tn k h
tn1
xm cnh
xm1 xm
xm+1
Supposons de plus que est obtenu de faon exprimentale et que lon mesure par exemple c e (xm + ph) = (xm + ph) + (1)p , i.e. on commet une erreur de sur || . n n Dans ce cas, lerreur sur vm est de lordre de (1 + 2c)n et, pour x, lerreur sur vm e cro de faon exponentielle avec le nombre ditrations n. On dit que le schma (5.1) est t c e e instable. On constate sur cet exemple que, bien quil semble facile de remplacer les drives partielles e e par des dirences nies, il faut toujours garantir que la solution donne par un schma e e e aux dirences nies converge, en un sens a dnir, vers une solution de lquation aux e ` e e drives partielles. e e
5.1. INTRODUCTION
69
Remarques : 1. Dans notre discussion du schma discret (5.1) on na pas tenu compte dventuelles e e conditions au bord, pour un schma discret celles-ci peuvent tre de deux types. e e Supposons que lon consid`re lquation de transport sur R+ R uniquement. e e + Pour pouvoir rsoudre lquation aux drives partielles on se donne une condition e e e e en x = 0 de la forme u(0, t) = (t) pour t 0 , ceci entra pour le schma discret ne e n que v0 = (tn ) (n 0) . De plus, tant donne la mmoire limit dun ordinateur, on est oblig de se rese e e e e treindre a une grille spatiale nie {x0 , x1 , . . . , xM } . En xM on introduit alors des ` n+1 n+1 n+1 n contraintes de la forme vM = vM 1 ou vM = vM 1 , . . . n+1 Ces conditions au bord numriques sont ncessaires pour pouvoir dterminer v m e e e pour 0 m M . Dans ce qui suit nous ninsisterons pas sur les probl`mes poss e e par ces contraintes. 2. Les schmas aux dirences nies permettent de traiter facilement des probl`mes en e e e une, deux ou trois dimensions spatiales pour des domaines a gomtrie simple et avec ` e e des proprits (p.ex. densit, vitesse, . . .) qui varient lentement dans le domaine de ee e dnition. D`s que la gomtrie devient complexe ou que les proprits changement e e e e ee rapidement, on prfrera la mthode des lments nis. ee e ee En eet, la mthode des lments nis permet de traiter des gomtries complexes e ee e e (e.g. bords internes) et la triangulation peut tre localement rane pour traiter e e des variations rapides de la solution. Par contre, il est en gnral plus compliqu e e e dcrire un code pour les lments nis que pour les schmas aux dirences nies. e ee e e
5.2
5.2.1
Dans toute la suite L est un oprateur direntiel linaire qui est dordre un en t et on e e e suppose que le probl`me e Lu(x, t) = f (x, t) dans R + u(x, 0) = (x) pour x ouvert de R (5.2)
est bien pos (voir page 3). e En remplaant les drives partielles par des dirences nies on obtient un oprateur c e e e e discret Lh,k . Ainsi lquation aux drives partielles homog`ne discrtise scrit sous la e e e e e e e forme Lh,k v = 0 . Pour tenir compte du second membre on utilise un oprateur discret e Ih,k que lon applique a f . ` Un schma aux dirences nies gnral, associ a lquation aux drives partielles (5.2), e e e e e` e e e scrit donc comme suit : e Lh,k v = Ih,k f .
n Dans la suite on choisira Ih,k toujours de faon a avoir Ih,k f = fm , de plus la condition c ` 0 initiale sera simplement donne par vm = (xm ) , pour xm . e Une premi`re classication des schmas discrets est donn par la dnition suivante. e e e e
70
Dfinition 5.2.1 e n+1 Un schma aux dirences nies est explicite si on peut crire vm comme combinaison e e e j linaire nie des vi pour j n . e n+1 e Le schma est dit implicite si dautres valeurs de v sont ncessaires (p.ex. vm1 ). e Un schma aux dirences nies est a un pas de temps sil utilise les valeurs de v a e e ` ` deux instants seulement, par exemple en tn et tn+1 . Le schma est a pas multiples si la valeur de v a plus de deux instants intervient. e ` `
n Un schma a un pas construit (vm )m , pour tout n 1 , a partir des conditions initiales e ` ` 0 0 vm = (xm ) . Pour initialiser un schma a J pas, les vm ne sont pas susants, il faut e ` j fournir (vm )m pour J instants.
Dans la section prcdente on a montr quil faut sassurer que la solution discr`te obtenue e e e e par un schma aux dirences nies donne une reprsentation correcte de la solution de e e e lquation aux drives partielles, do` la e e e u Dfinition 5.2.2 (Convergence) e n Soit u(x, t) la solution de (5.2) et v une solution du schma discret Lh,k v = fm telle que e 0 vm converge vers (x) quand xm converge vers x . n On dit que le schma aux dirences nies Lh,k est un schma convergent si vm converge e e e vers u(x, t) quand (xm , tn ) converge vers (x, t) pour (h, k) tendant vers (0, 0) . La convergence exprime que la solution dun schma aux dirences nies converge vers la e e solution de lquation aux drives partielles. Avant de lancer des calculs, il faut sassurer e e e que le schma utilis est convergent, or ceci est en gnral dicile a dmontrer. On va e e e e ` e proposer dans la suite des crit`res plus facile a vrier. e ` e Remarques : 1. La dnition ci-dessus et celles qui suivent, ne concernent que les schmas a un e e ` pas de temps. Pour des mthodes a pas multiples il faudra tenir compte de ltape e ` e dinitialisation.
n 2. La dnition prcise de la convergence de vm vers u , cohrente avec la suite, entra e e e ne n 2 lutilisation dun oprateur dinterpolation : pour n x, on associe a (vm )m l (hZ) e e ` un lment de L2 (R) . On impose ensuite une convergence dans L2 (R) de linterpole ee e vers la solution u .
Dfinition 5.2.3 (Consistance) e n e e e On dit que le schma Lh,k v = fm est consistant avec lquation aux drives partielles e Lu = f si pour toute fonction de classe C on a, en tout point (xm , tn ) :
(h,k)(0,0)
lim
(L Lh,k ) = 0 .
La consistance entra en particulier quune solution rguli`re de lquation aux drives ne e e e e e partielles est une solution du schma aux dirences nies quand les pas de discrtisae e e tion tendent vers 0 . Cest une condition ncessaire de convergence mais ce nest pas une e condition susante. En eet, considrons de nouveau le schma (5.1) et soit une fonction de classe C . e e n+1 n n n m m + c m+1 . On a L = t + cx et Lh,k = m k h
71
Grce a la formule de Taylor : a ` n m+1 n+1 m h2 = (xm + h, tn ) = + hx (xm , tn ) + xx (xm , tn ) + O(h3 ) , 2 k2 = (xm , tn + k) = n + kt (xm , tn ) + tt (xm , tn ) + O(k 3 ) . m 2 n m
On en dduit que e 1 Lh,k (xm , tn ) = t (xm , tn )+c x (xm , tn )+ (ktt (xm , tn ) + c hxx (xm , tn ))+O(h2 )+O(k 2 ) 2 et donc 1 (L Lh,k )(xm , tn ) = (ktt (xm , tn ) + c hxx (xm , tn )) + O(h2 ) + O(k 2 ) (h,k)(0,0) 0 . 2 Le schma (5.1) est donc consistant bien que non convergent. e On avait constat que ce schma est en un certain sens instable. Nous allons donner une e e dnition prcise de stabilit pour une quation aux drives partielles homog`ne dordre e e e e e e e un en t dans la Dfinition 5.2.4 (Stabilit) e e Le schma aux dirences nies Lh,k v = 0 , associ a lquation aux drives partielles e e e ` e e e 2 Lu = 0 , est stable sil existe (R+ ) , avec (0, 0) , tel que pour tout T > 0 il existe une constante CT tel que
+ + n |vm |2
h
m=
CT h
m=
0 |vm |2 ,
pour tous 0 tn T et (h, k) . On dit que est la rgion de stabilit du schma. e e e Un exemple frquent pour est le segment {(h, h) / 0 < h < C} , avec C et des e constantes positives ; le cne {(h, k) / 0 < C1 h k C2 h < C3 } , avec C1 , C2 et C3 des o constantes positives est un autre exemple. Remarque : En utilisant le principe de Duhamel (voir p. 42) on montre quun schma e n Lh,k v = fm est stable si le schma homog`ne Lh,k v = 0 est stable. De plus on peut tendre e e e la dnition au cas dquations dordre 2 en t et au cas de schmas a pas multiples. e e e ` La dnition de stabilit utilise une norme sur lespace l 2 (hZ) (voir section 5.3) : e e
+ n (vm )mZ 1/2 n |vm |2
l2 (hZ)
h
m=
est ni.
Le crit`re de stabilit scrit alors, pour h et k susamment petits dans : e e e (T > 0)(CT )(tn [0, T ]) : vn
l2 (hZ)
CT v 0
l2 (hZ)
(S)
La stabilit garantit qu` chaque instant tn [0, T ] , la norme de la solution discr`te est e a e borne, a un facteur constant pr`s, par la norme des donnes initiales. e ` e e Cette stabilit numrique nest pas a confondre avec la stabilit du probl`me bien pos e e ` e e e (5.2). Cette derni`re concerne le comportement a linni de la solution en fonction des e `
72
conditions initiales, par contre la stabilit numrique dun schma concerne le comportee e e ment de v sur lintervalle [0, T ] quand (h, k) tend vers (0, 0) . Remarques : 1. Il existe dautres faons de dnir la stabilit numrique, on a choisi la dnition c e e e e particuli`re ci-dessus parce quelle sapplique a un grand nombre de schmas linaires e ` e e et quelle est adapte a lanalyse de stabilit de von Neumann . e ` e 2. Pour des quations aux drives partielles et schmas discrets non linaires on utilie e e e e 1 sera dautres espaces, e.g. L , et dautres proprits, e.g. la diminution de la variaee tion totale ou la monotonie dun schma, an dviter des oscillations de v l` o` la e e a u solution de lquation aux drives partielles varie rapidement, e.g. discontinuits. e e e e Limportance des notions de consistance et stabilit est due au e Thor`me 5.2.1 (Equivalence de Lax-Richtmyer) e e Un schma linaire, consistant avec le probl`me (5.2) , est convergent si et seulement si il e e e est stable. La derni`re dnition que lon va introduire permet de distinguer les schmas convergents, e e e en eet deux schmas convergents nont pas ncessairement les mmes performances pour e e e lapproximation de la solution continue. Dfinition 5.2.5 e n Un schma Lh,k v = fm , consistant avec le probl`me (5.2) , est dordre p en espace et q en e e temps si, pour toute fonction de classe C : Eh,k = Lh,k L = O(hp ) + O(k q ) . On dit que le schma est dordre (p, q) et Eh,k est lerreur de troncature du schma. e e Remarques : 1. Lerreur de troncature permet dvaluer lerreur introduite lorsque lon remplace e loprateur continu L par loprateur discret Lh,k . e e 2. La consistance dun schma ncessite uniquement que Lh,k L = o(1) . e e Exemple : Sans dmontrer le thor`me de Lax-Richtmyer nous allons illustrer, grce a un exemple, e e e a ` comment la stabilit et la consistance interviennent dans la convergence dun schma e e discret. Pour c < 0 , considrons lquation dadvection ut + cux = 0 dans R R , e e + avec u(x, 0) = (x) pour x R . On suppose de plus que est borne sur R . e n+1 n n On va utiliser le schma (5.1) : vm = (1 + c)vm cvm+1 , o` = k/h . e u n Pour x = mh, posons e(x, tn ) = u(x, tn ) vm , cest lerreur commise par le schma au e point (x, tn ) . Alors Lh,k e = Lh,k u = Eh,k u , car v est solution de Lh,k v = 0 et u est solution de Lu = 0 , do` u e(x, tn+1 ) = (1 + c) e(x, tn ) c e(x + h, tn ) + k Eh,k u(x, tn ) et sup |e(x, tn+1 )| (|1 + c| + |c|) sup |e(x, tn )| + k sup |Eh,k u(x, tn )| .
x x x
73
On a ainsi major lerreur a linstant tn+1 grce a lerreur a linstant tn et de lerreur de e ` a ` ` troncature en tn ; par rcurrence e
n
p=1
Or comme lon verra plus loin, lanalyse de von Neumann permet de montrer que le schma est stable si et seulement si est tel que 1 c 0 , i.e. |1 + c| + |c| = 1 . e Comme de plus e(x, 0) = 0 pour tout x , on a
n1
sup |e(x, tn )| k
x
i=0
On doit maintenant valuer lerreur de troncature. Or le schma (5.1) est consistant et e e dordre (1, 1) , de plus Eh,k u(x, tn ) = Lh,k u(x, tn ) = (Lh,k uLu)(x, tn ) = O(k) car k = h . On obtient ainsi une estimation de lerreur de troncature qui dpend de uxx , utt et k , e pour obtenir une estimation globale et montrer la convergence du schma on va utiliser e lhypoth`se sur . On a e |Eh,k u(x, ti )| = |Lh,k (x cti )| = et nalement, pour tout n T /k : e(., tn )
k h c (1 ) + c2 (2 ) k C sup | ()| 2 2
,c)
T C(,
,c)
k = C(,
,c,T ) k
Lerreur commise en chaque point de la grille est dordre k et le schma est convergent. e Remarque : Dans cet exemple nous avons montr une convergence en chaque point de e la grille en utilisant la norme . sur hZ et des hypoth`ses fortes sur u(x, 0) . e Nous allons appliquer lquivalence de Lax-Richtmyer uniquement dans le cadre L 2 , en e eet on a alors, grce a la mthode de von Neumann, une mthode gnrale pour ltude a ` e e e e e de la stabilit pour des schmas numriques associs a des quations aux drives pare e e e ` e e e tielles dvolution linaires a coecients constants. e e ` On montre par ailleurs que la rgularit des donnes initiales a une inuence sur la pre e e e cision de la solution numrique. e
5.2.2
Condition de Courant-Friedrichs-Lewy
Nous allons tudier dirents schmas discrets pour le probl`me de transport avec la e e e e vitesse c R : ut + c u x = 0 dans R R + u(x, 0) = (x) pour x R (5.4)
n+1 n n vm = vm + vm+1
=h
mZ
n+1 |vm |2 = h
mZ
2 l2 (hZ)
74
Une condition ncessaire de stabilit de ces schmas est donc || + || 1 . e e e En particulier pour le schma (5.1) on a |1+c|+|c| 1 , ce schma est donc stable pour e e 1 c 0 . On va montrer plus loin que ceci est une condition susante de stabilit. e Le schma de Lax-Friedrichs est centr en espace et progressif en temps mais utilise e e n une moyenne centre pour approximer vm : e
1 n n n+1 n vm 2 (vm+1 + vm1 ) v n vm1 + c m+1 =0. k 2h
(5.5)
n+1 n n Cest un schma consistant et qui est de la forme vm = vm+1 + vm1 . e On montre que ces schmas sont stables si || + || 1 , en particulier le schma de e e Lax-Friedrichs (5.5) est stable si |c| 1 .
Les rsultats prcdents se gnralisent et on a le e e e e e Thor`me 5.2.2 e e n+1 n n n Soit vm = vm1 + vm + vm+1 un schma explicite pour lquation aux drives e e e e partielles (5.4) . Alors, si le rapport k/h = est constant, une condition ncessaire de e stabilit du schma est la condition de Courant-Friedrichs-Lewy (CF L) : e e Thor`me 5.2.3 e e Il nexiste pas de schma explicite, consistant et inconditionnellement stable pour lquae e tion aux drives partielles (5.4) . e e Remarque : 1. Les thor`mes prcdents sappliquent de faon plus gnrale aux syst`mes dquae e e e c e e e e tions aux drives partielles hyperboliques, dnis a la page 45 . e e e ` 2. La vitesse dadvection c dtermine la pente des droites caractristiques de (5.4), la e e condition (CF L) assure que le domaine de dpendance discret contient le domaine e de dpendance de lquation aux drives partielles, voir la gure page 68. e e e e Exemple : Pour montrer que le dernier thor`me ne sapplique pas aux schmas implicites e e e considrons le schma suivant, progressif en temps et rtrograde en espace, a linstant t n+1 : e e e `
n+1 n+1 n v n+1 vm1 vm v m +c m =0. k h n+1 n n+1 En crivant ce schma sous la forme : (1 + c)vm = vm + cvm1 e e on montre que si c > 0 , on a pour tout R+ :
|c| 1 .
v.n+1
2 l2 (hZ)
v.n
2 l2 (hZ)
5.3
Lanalyse de stabilit de von Neumann utilise lanalyse de Fourier pour tudier le come e portement dun schma discret. Avant de montrer des applications a ltude de schmas e ` e e discrets, nous allons motiver lutilisation de lanalyse de Fourier grce a une analogie et a ` 2 nous allons rappeler quelques rsultats sur l (hZ) . e
75
Motivations : Si on consid`re que v.n est un signal numrique discret monodimensionnel on peut intere e prter v.n+1 comme la sortie dun syst`me discret S : v.n+1 = S[v.n ] . e e Pour tudier le syst`me S on introduit la transforme de Fourier des signaux discrets v n e e e et v n+1 . Ceci est particuli`rement intressant quand S est un syst`me linaire et invariant, un e e e e ltre, dans ce cas la sortie de S est obtenue grce a la convolution de lentre du syst`me a ` e e avec la rponse impulsionnelle s de S : S[e] = s e . Un ltre peut aussi tre reprsent e e e e par une quation aux dirences nies, rcursive ou non. e e e Dans le domaine frquentiel la convolution devient une multiplication de la transforme e e de Fourier du signal en entre avec la rponse frquentielle de S . Le signal sinuso e e e s dal im discret e = (e )mZ est un vecteur propre de S : S[e ] = 2 s() e . Lespace l2 (hZ) : Nous rappelons ici bri`vement quelques lments de la structure de lespace de Hilbert e ee l2 (hZ) . Les dnitions et rsultats prsents sobtiennent a partir de ltude des sries e e e e ` e e 2 de Fourier des fonctions de L ([, +]) : grce a un changement de variable, on tient a ` compte du pas de discrtisation spatiale h > 0 . e
+ 1/2
l2 (hZ)
h
m=
|vm |
vm wm .
Soit v l2 (hZ) , on dnit sa transforme de Fourier, pour [/h, +/h] , par : e e 1 v() = 2
+
hvm eimh .
m=
La fonction v est 2/hpriodique, cest un lment de L2 ([/h, +/h]) et lidentit e ee e ci-dessus est a prendre au sens de la convergence en moyenne quadratique. ` La formule de reconstruction scrit pour tout m Z : e 1 vm = 2 et on a lidentit de Parseval : e v
2 L2 ([/h,/h])
+h
v()e+imh d ,
+ h h
|v()| d =
m=
h|vm |2 = v
2 l2 (hZ)
n n Pour v l2 (hZ) on a : Th v() = e ih v() , o` Th vm = vm 1 ; u une variation spatiale (translation) de v se traduit par un dphasage de v . e
76
Grce a la formule de Parseval on obtient une forme quivalente pour le crit`re de a ` e e stabilit (S) : e (T > 0)(CT )(tn [0, T ]) : Mthode de von Neumann : e Considrons un schma discret a un pas Lh,k v = 0 , linaire et a coecients constants. En e e ` e ` appliquant la transforme de Fourier et en rarrangeant les termes on obtient la relation : e e v n+1 () = g(, h, k) v n() . Le facteur g(, h, k) est appel facteur damplication : le module de g , |g| , est lamplie cation et la phase de g , arg(g) , est le dphasage que subit chaque frquence de la solution e e discr`te quand on avance dun pas de temps k . Par itration on a : e e v n () = g(, h, k)n v 0 () . Remarques : 1. Pour dterminer le facteur damplication en pratique, on remplace dans le schma e e n n imh n imh discret les termes vm par g e , i.e. on cherche des solutions de la forme g e . 2. Si u vrie (5.4) et, pour tout t R+ , u(., t) L2 (R) , on a, en appliquant la e transforme de Fourier par rapport a x (voir aussi la section 2.3.2) : e ` ut (, t) = c i u(, t) , donc u(, t) = u(, 0) eict = () eict et u(, t + k) = eick u(, t) . Le facteur damplication, g = |g| ei arg(g) , du schma, doit approximer eick : e Si |g| = 1 on a une erreur damplitude dans la solution discr`te, en particulier si e |g(, h, k)| < 1 , on attnue la frquence et on parle de dissipation, resp. viscosit, e e e numrique ; la dirence de phase, arg(g) ck , introduit une erreur de phase, e e appele dispersion numrique . e e Exemples : 1. Considrons lquation de transport ut + c ux = 0 , pour c > 0 , discrtise grce au e e e e a schma progressif en temps et en espace (5.1) . e On a g(, h, k) = 1 + c c eih qui, pour x, ne dpend que de h , et e e |g(h)|2 = g(h) g(h) = 1 + 4c(1 + c) sin2 (h/2) . Or |v n ()| = |g(h)|n |v 0 ()| , v.n
L2 ([/h,/h])
CT v.0
L2 ([/h,/h]))
(S)
et |g(h)| = 1 seulement si h = 0 ; sinon |g(h)| > 1 , on va en dduire que le e schma est instable. e Soit (h0 , k0 ) (R )2 , et K R : + + alors il existe (h, k) ]0, h0 ]]0, k0 ] et (1 , 2 ) ]0, /h[2 avec |g(h)| 1 + Kk , pour tout [1 , 2 ] .
5.3. ANALYSE DE STABILITE DE VON NEUMANN
77 v0
l2 (hZ)
(2 1 )1/2 si [1 , 2 ] et 0 si [1 , 2 ]
= 1.
= =
vn
2 L2 ([/h,/h])
+ h h 2
=
1
|g(h)|2n
(1 + Kk)2n
pour n T /k et k0 susamment petit. Comme K peut tre pris aussi grand que lon veut, le crit`re de stabilit (S) nest e e e pas vri et le schma (5.1) est instable, donc divergent, pour c > 0 . e e e 2. Considrons lquation aux drives partielles ut + c ux u = 0 . e e e e Ce probl`me dadvection-raction est discrtis grce a (5.5) : e e e e a `
1 n n n+1 n vm 2 (vm+1 + vm1 ) v n vm1 n + c m+1 vm = 0 . k 2h
Le facteur damplication est g(, h, k) = cos(h) + i c sin(h) + k , et |g|2 = (cos(h) + k)2 + c2 2 sin2 (h) . On montre que |g(, h, k)| 1 + k pour |c| 1 , or pour n T /k : vn
2 l2 (hZ)
= vn
2 L2 ([/h,/h]) 2 L2 ([/h,/h])
(1 +
(1 + k)2n v 0
= |g(, h, k)|2n v 0
e2T v 0
2 l2 (hZ)
Dapr`s le crit`re (S) on a la stabilit du schma. e e e e La solution de lquation aux drives partielles scrit u(x, t) = (x ct) et , e e e e ainsi mme si la condition initiale est borne sur R , la solution u nest pas e e borne sur R R+ . En appliquant la transforme de Fourier en espace, on a e e u(, t) = () eict et et donc u(, t + k) = u(, t) eick ek . An de pouvoir suivre lvolution de la solution u(x, t) au cours du temps, il faut e que |g| soit plus grand que 1.
78
De faon gnrale on a le c e e Thor`me 5.3.1 e e Un schma a un pas Lh,k v = 0, linaire et a coecients constants, est stable si et seulement e ` e ` si il existe une constante K R+ et (R )2 , avec (0, 0) , tel que + pour tout (h, k) . Si le facteur damplication ne dpend que de h la condition ncessaire et susante de e e stabilit scrit : e e |g(h)| 1 . Exemple de Strikwerda : On consid`re lquation aux drives partielles ut + c uxxx = f , e e e e o` lon discrtise ut comme dans le schma de Lax-Friedrichs (5.5) et uxxx grce aux u e e a dirences centres, pour obtenir : e e
1 n n n+1 n n n n vm 2 (vm+1 + vm1 ) vm+2 2vm+1 + 2vm1 vm2 n +c = fm . 3 k 2h Soit une fonction rguli`re, alors L = t + c xxx et e e
|g(, h, k)| 1 + Kk ,
n+1 1 (n + n ) n 2n + 2n n m m+1 m1 m+1 m1 m2 2 + c m+2 . k 2h3 En utilisant la formule de Taylor en (xm , tn ) : Lh,k = n+1 = n + k t + m m n m1 n m2 Donc Lh,k = t + et k2 tt + O(k 3 ) , 2 h3 h4 h2 xx xxx + xxxx + O(h5 ) , = n h x + m 2 6 4! 4h3 2h4 = n 2h x + 2h2 xx xxx + xxxx + O(h5 ) . m 3 3 h2 k tt xx + O(k 2 + h4 /k) + cxxx + O(h2 ) 2 k
k h2 tt xx + O(k 2 + h4 /k + h2 ) . 2 k Le schma est donc consistant seulement si h2 /k tend vers 0 quand h, k 0 . e Supposons que k = H(h) , avec H rguli`re et H(0) = 0 , dans ce cas Lh,k L = O(h) e e et on dit que le schma est dordre 1 . e Lh,k L = Pour tudier la stabilit on crit le schma sous la forme : e e e e ck n 1 n n n n n n+1 vm = (vm+1 + vm1 ) 3 (vm+2 2vm+1 + 2vm1 vm2 ) , 2 2h on obtient le facteur damplication g(, h, k) = cos(h) + i 4ckh3 sin(h) sin2 (h/2) . On a |g(, h, k)|2 = cos2 (h) + 16c2 k 2 h6 sin2 (h) sin4 (h/2) , une condition ncessaire de stabilit est alors que la quantit h3 k reste born. e e e e 2 1 Or ceci est incompatible avec la contrainte de consistance sur h k . Le schma ne peut e donc pas tre consistant et stable en mme temps. e e
79
Schmas pour lquation de transport. e e Pour terminer on prsente quelques schmas classiques pour la rsolution du probl`me e e e e homog`ne (5.4) pour c R : e Schma centr en espace : e e
n n n+1 n vm+1 vm1 vm v m +c =0; Ce schma scrit e e k 2h il est consistant et dordre (2, 1) . En crivant e n n+1 vm = v m
c n n vm1 ) , (v 2 m+1
o` = k/h , on montre que le facteur damplication g(h) = 1 i c sin(h) u et |g|2 = 1 + c2 2 sin2 (h) 1 : on a un schma instable. e
Schma e
upwind
On a vu que pour c > 0 il faut utiliser une dirence nie rtrograde pour garantir la e e stabilit et pour c < 0 il faut utiliser une dirence nie progressive. En intgrant le test e e e n n sur le signe de c , on obtient un schma qui utilise vm1 ou vm+1 , suivant la direction e dadvection, do` le nom du schma. On a u e
n+1 n vm v m + k
c |c| 2
n n vm+1 vm + h
c + |c| 2
n n vm vm1 = 0; h
On obtient donc un schma explicite centr en espace avec un terme correcteur provenant e e dune discrtisation de (|c|h/2) uxx . e Le schma est consistant, dordre (1, 1) et g(h) = 1 2|c| sin2 (h/2) i c sin(h) . e On a |g(h)|2 = 1 4|c|(1 |c|) sin2 (h) 1 si et seulement si |c| 1 . On constate que le terme de dissipation ou viscosit numrique (|c|h/2) u xx a rendu stable e e le schma centr en espace, par contre on a perdu dans lordre du schma. e e e Schma de Lax-Friedrichs : e Ce schma se prsente sous la forme suivante e e
n n+1 n n vm 1 (vm+1 + vm1 ) v n vm1 2 + c m+1 =0. k 2h
Pour rguli`re on montre que e e Lh,k L = h2 ch2 k tt xx + O(k 2 + h4 /k) + xxx + O(h4 ) . 2 2k 6
80
Si k = h , avec R x, on a un schma consistant dordre 1. e e + Le facteur damplication est g(h) = cos(h) + i c sin(h) et |g|2 = cos2 (h) + c2 2 sin2 (h) 1 si et seulement si |c| 1 . On peut crire ce schma sous la forme e e
n n n n n n+1 n vm+1 vm1 1 vm+1 2vm + vm1 vm v m +c =0 k 2h 2 k
et on constate que le schma de Lax-Friedrichs correspond a la discrtisation par die ` e e rences nies centres de lquation aux drives partielles e e e e ut + c u x h2 uxx = 0 . 2k
On observe encore le phnom`ne de viscosit numrique qui rend le schma stable. e e e e e Schma de Lax-Wendroff : e Ce schma est bas sur lutilisation du dveloppement de Taylor et de lquation aux e e e e drives partielles a discrtiser. Si u est une solution rguli`re de lquation aux drives e e ` e e e e e e partielles on a : u k 2 2 u + + O(k 2 ) , un+1 = un + k m m t 2 t2 or ut = cux et utt = ( cux )t = cuxt = c2 uxx , do` u un+1 = un ck m m k2 2 u u + c2 + O(k 2 ) . x 2 x2
En utilisant une discrtisation spatiale centre pour approximer ux et uxx , on obtient le e e schma discret : e
n+1 vm
n vm
Ce schma est consistent et dordre (2, 1) car e c2 k c2 kh2 k xx + xxxx + O(h3 + k 2 ) . Lh,k L = tt 2 2 4! En posant = k/h :
n+1 n vm = v m
on montre que le facteur damplication est g(h) = 1 2c2 2 sin2 (h/2) i c sin(h) et |g|2 = 1 4c2 2 (1 c2 2 ) sin4 (h/2) . Ce schma est donc stable pour la condition (CF L) : |c| 1 . e
81
ou
n+1 n1 n n vm = vm c (vm+1 vm1 ) . 0 1 Ce schma est a deux pas, il faut donc donner les conditions initiales en vm et vm . En e ` pratique on utilise un schma a un pas pour initialiser un schma a pas multiples. e ` e ` 2 2 k ch Comme Lh,k L = ttt xxx + O(h2 + k 2 ) , 6 6 le schma est consistant et dordre (2, 2) . e
Lanalyse de von Neumann ne sapplique pas directement a ce schma. En appliquant la ` e transform de Fourier spatiale on a e v n+1 () + i2c sin(h) v n () v n1 () = 0 . (5.7)
Cette rcurrence dordre deux sexprime grce a un syst`me en introduisant une variable e a ` e supplmentaire : e v n+1 () v n () Si on pose V n = vn v n1 = i2c sin(h) 1 1 0 v n () v n1 ()
Si le rayon spectral (G) = max{|gi | / gi valeur propre de G} est strictement plus petit que 1, alors on sait que Gn tend vers 0 et Gn est borne. e Mais si lon veut gnraliser le cas scalaire, on sait quil faut permettre |gi | 1 + Kk . Or e e pour les syst`mes ceci est une condition ncessaire de stabilit mais non susante, comme e e e on le verra dans ce qui suit. Nous allons tudier les valeurs propres de G . Remarquons dabord que G est la matrice e compagne de la relation de rcurrence (5.7), son polynme caractristique correspond a e o e ` n lquation en g, obtenue en remplaant vm par g n eimh dans le schma discret : e c e g 2 + i2c sin(h) g 1 = 0 . e Cas c < 1 : on a deux racines complexes conjugues g+ et g avec |g | = 1 . La matrice damplication G(h) est diagonalisable et la solution de (5.7) est v n () = (h) g+ (h)n + (h) g (h)n , o` et dpendent des donnes initiales v 0 et v 1 . On montre que u e e vn =
n n n n g+ g g + g 0 g+ g 1 g+ v 0 v 1 g v 0 + v 1 g+ + g = v + v . g+ g g+ g g+ g g+ g
Schma e
leapfrog
(saute-mouton) :
G n CT .
1 c2 2 sin2 (h) .
82 On en dduit que e
2 1 c 2 2
v0
l2 (hZ)
+ v1
l2 (hZ)
Le schma est donc stable : a tout instant tn [0, T ] , la norme de v n est borne par celle e ` e des conditions initiales. Cas |c| > 1 : si h = /2 on a deux racines distinctes g (/2) . Or |g+ (/2)|2 = 1 + 2c2 2 > 1 : la matrice G(/2) est diagonalisable mais Gn (/2) nest pas born et le schma nest pas stable. e e Cas |c| = 1 : on a des racines doubles pour h = /2 et la solution de (5.7) scrit e v n () = (h) g+(h)n + (h) ng+ (h)n1 , o` et dpendent de v 0 et v 1 . On a u e v n () = g+ (h)n (1 n) v 0 () + ng+ (h)n1 v 1 () . Or, pour h = /2 , on trouve g+ (h) = i et, si lon suppose que v 1 () = g0 (h) v 0 () , on a v n = (i)n1 (ng0 i(1 n)) v 0 .
Pour = 2h on a donc |v n ()| nC |v 0 ()| pour n grand, on peut donc construire des solutions pour lesquelles v n l2 (hZ) cro de faon linaire avec n . t c e
Dans ce cas la matrice damplication G(/2) a une valeur propre double de module 1 et elle peut se mettre sous forme de Jordan mais le schma est instable dapr`s la thorie. e e e On peut remarquer que cette instabilit est tr`s faible et rarement nuisible en pratique. e e
n+1 n+1 Note : On remarque que les valeurs de vm et vm+1 sont indpendantes, on pourra voir e appara deux solutions indpendantes (en chiquier), ceci peut poser des probl`mes en tre e e e pratique mme si |c| < 1 (en particulier pour des probl`mes non linaires). e e e
83
5.4
Equations paraboliques
Les dnitions de convergence, consistance, stabilit et dordre dun schma sappliquent e e e aussi aux quations paraboliques. Notre quation type sera e e ut = d uxx , Schma explicite centr en espace : e e On a
n n n n+1 n vm1 2vm + vm+1 vm v m d = 0, k h2
avec d R . +
et Lh,k L =
k dh2 tt + O(k 2 ) xxxx + O(h4 ) , 2 12 le schma est consistant et dordre (2, 1) . En posant = k/h2 : e
n+1 n n n vm = d (vm1 + vm+1 ) + (1 2d) vm
et le facteur damplication g(h) = 1 4d sin2 (h/2) . Donc |g| 1 si et seulement si 0 4d sin2 (h/2) 2 pour tout . La condition de stabilit est donc e d 1 . 2
On constate que la condition de stabilit entra que k/h 0 quand h, k 0 . Ceci e ne exprime le fait que la solution discr`te v doit avoir un domaine de dpendance proche de e e celui de la solution u de lquation aux drives partielles qui est ici laxe {t = 0} . e e e Lquation de la chaleur modlisant une diusion il est intressant davoir un schma e e e e dissipatif, i.e. |g| < 1 on choisira donc < 1/(2d) . 1 Notons nalement que pour ]0, 2d [ x et k = h2 le schma est dordre 2. e e Schma implicite centr en espace : e e En crivant un schma rtrograde en temps et centr en espace on a : e e e e
n+1 n+1 n+1 n v n+1 2vm + vm+1 vm v m d m1 = 0, k h2
ou,
n+1 n+1 n+1 n (1 + 2d) vm d (vm1 + vm+1 ) = vm .
` A chaque tape il faut rsoudre un syst`me tridiagonal pour obtenir v n+1 a partir de v n . e e e ` Ce schma est consistant, dordre (2, 1) et inconditionnellement stable car e g(h) = 1 . 1 + 4d sin2 (h/2)
84
Schma de Crank-Nicolson e Ce schma est base sur les deux prcdents, on value le terme de diusion en faisant la e e e e e moyenne de lcriture explicite et implicite : e
n+1 n vm v m d k 2 n+1 n+1 n n n n+1 vm1 2vm + vm+1 vm1 2vm + vm+1 + h2 h2
= 0.
1 2d sin2 (h/2) 4d sin2 (h/2) =1 , 1 + 2d sin2 (h/2) 1 + 2d sin2 (h/2) donc |g| 1 et le schma est inconditionnellement stable. e Le facteur damplication est g(h) = Schma de Du Fort-Frankel e Le schma leapfrog est instable pour tout , or une modication de la discrtisation e e spatiale donne le schma stable suivant e
n+1 n1 n n+1 n+1 v n (vm + vm ) + vm+1 vm v m d m1 = 0, 2k h2
avec Lh,k L = dk 2 h2 tt + O(h2 ) + O(k 2 ) . Ce schma est donc consistant seulement si k/h 0 quand h, k 0 . En notant e = k/h2 on a
n+1 n n n1 (1 + 2d) vm = 2d (vm+1 + vm1 ) + (1 2d) vm .
Ltude de la stabilit nous am`ne a considrer le polynme caractristique de la matrice e e e ` e o e damplication G : (1 + 2d) g 2 4d cos(h) g (1 2d) = 0 . Dont les racines scrivent : g (h) = e 2d cos(h) 1 4d2 2 sin2 (h) . 1 + 2d
Si 4d2 2 sin2 (h) < 1 , on a deux racines relles distinctes et e |g (h)| 2d| cos(h)| + 1 4d2 2 sin2 (h) <1. 1 + 2d
Si 4d2 2 sin2 (h) > 1 , on a deux racines complexes conjugues et e |g (h)|2 = 1 4d2 2 1 (2d cos(h))2 + 4d2 2 sin2 (h) 1 = <1. (1 + 2d)2 (1 + 2d)2 |g+ (h)| =
2d| cos(h)| <1. 1 + 2d Dans tous les cas le rayon spectral est strictement plus petit que 1 et G(h)n tend vers zro. e On a donc un schma explicite qui est inconditionnellement stable mais conditionnellement e consistant.
85
5.5
5.5.1
Applications
Equation dadvection-diusion-raction e
On sintresse a la rsolution numrique de lquation aux drives partielles suivante : e ` e e e e e u u 2u (x, t) + c (x, t) = d 2 (x, t) + r u(x, t) , t x x pour (x, t) R R + (5.8)
(5.9)
1 2
Ce schma est consistant et dordre (2, 1) , la condition de stabilit est d e e 1 e e Pour ]0, 2d [ x et k = h2 on a un schma dordre 2. ch , le schma (5.9) scrit aussi e e En notant = 2d
n+1 n n n vm = (1 2d + rk)vm + d(1 )vm+1 + d(1 + )vm1 . 1 Pour ]0, 2d [ x, 1 et 1 2d + rh2 0 , on a e
Pour obtenir ce rsultat, on a impos des contraintes sur la discrtisation spatiale h : e e e h 2d c et, si r < 0, k = h2 2d 1 . r
Si ces ingalits ne sont pas vries, la solution discr`te na pas le mme comportement e e e e e e que la solution de lquation aux drives partielles parabolique (5.8) et on pourra observer e e e des oscillations dans la solution. Par contre, si pour satisfaire ces contraintes il faut prendre h trop petit, il se peut que le schma devient inecace. e Remarque : En modlisation dun ux de chaleur, la quantit Pe = c/d [1/m] est appele e e e nombre de Peclet , en mcanique des uides Re = c/d est le nombre de Reynolds. Ce e nombre contrle lintensit relative de ladvection par rapport a la diusion : si Pe est o e ` grand, cest ladvection qui domine, si Pe est petit cest la diusion. On peut interprter e 2 c/d comme tant le rapport du temps de diusion, tdif f = x /d , et du temps dadvection, e tadv = x/c . Ci-dessous on propose un code Scilab 1 pour implmenter le schma (5.9). Noter que lon e e sest intress principalement aux calculs et non pas a lutilisation du code (p.ex. entre e e ` e interactive des param`tres). e
1
86
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
// Param`tres de ledp e c=0.30; d=0.10; r=0.13; // Param`tres du schma e e k=0.1; T=10; N=ceil(T/k)+1; t=linspace(0,T,N); h=0.2; xg=-5; xd=10; M=ceil((xd-xg)./h)+1; x=linspace(xg,xd,M); v=zeros(N,M);
// // // // // // // //
pas de temps instant final nb. de points dans [0,T] intervalle discret de temps pas de discrtisation spatiale e limites spatiales nb. de points dans [xg,xd] intervalle discret despace
// Afficher les param`tres du schma e e mu=k/h^2; alpha=c.*h./(2.*d); mprintf(c=%f\td=%f\tr=%f\n,c,d,r) mprintf(h=%f\tk=%f\nmu=%f\talpha=%f \n,h,k,mu,alpha) mprintf(mu*d=%f\t2/Pe=%f\t,mu.*d,(2.*d)./c) if (r<0) then mprintf((2d*mu-1)/r=%f \n,(2.*d.*mu- 1)./r) else mprintf(\n) end // Condition initiale phi() dfinie sur [xg0,xd0] e xg0=-2;xd0=2; Ind_xg0_xd0=bool2s((xg0 <= x)&(x <= xd0)); // indicatrice de [xg0,xd0] deff([y]=phi(z),y=sin(%pi *z/2)); v(1,:)=Ind_xg0_xd0 .* phi(x); // dfinition de la c.i. e // Itrations sur le vecteur v(n,:) e c1=1-2.*d.*mu+r.*k; c2=d.*mu.*(1-alpha); c3=d.*mu.*(1+alpha); for n=2:N v(n,:)=c1*v(n-1,:) + c2*[v(n-1,2:M) 0] + c3*[0 v(n-1,1:M-1)]; end // Afficher la solution v dans [xg,xd]x[0,T]x[min(v),max(v)] plot3d1(x,t,v,-111,86,"x@t@u",[2 2 4])
On obtient lachage des param`tres : e c=0.300000 h=0.200000 mu=2.500000 mu*d=0.250000 d=0.100000 k=0.100000 alpha=0.300000 2/Pe=0.666667 r=0.130000
et, apr`s calculs, la gure suivante (on sest restreint a x [5, +5]) : e `
5.5. APPLICATIONS
87
u 1 t x 0 PSfrag replacements
5 1 10 5 0 5 0
On constate le transport et la diusion, a laquelle soppose la raction, de la concentration ` e initiale . Souvent on trace la solutions dans le plan xu a des instants tn dirents : ` e
1.0
t0
0.8 0.6 0.4 0.2
tn = nk avec k = 0.1
t25
t50
t75
t100
u
0.2 0.4
PSfrag replacements
0.6 0.8 1.0 5 4 3 2 1 0 1 2 3 4 5
Finalement, on sintresse a la prcision de v, solution de (5.9) sur la grille discr`te, par e ` e e rapport a u, solution de (5.8) sur R [0, T ] . On a a linstant t100 = 10, avec h = 0, 2 : ` ` err (t100 ) = et sup
xm [5,10] 100 u(xm , 10) vm = 0, 0497294
Sur la gure suivante on a reprsent la solution continue u et la solution discr`te v , on e e e 2 constate des dirences de lordre de 5 10 , cf. err . Lerreur relative, |u v|/|v|, peut e atteindre 2.6 aux points o` u 0 . u
CHAPITRE 5. SCHEMAS AUX DIFFERENCES FINIES
err2 (t100 ) = h
xm [5,10]
1/2
= 0, 0722882 .
88
Pour conclure on peut dire que le schma reproduit assez bien le comportement de la e solution de lquation aux drives partielles, par contre la discrtisation grossi`re et des e e e e e phnom`nes de bord perturbent rapidement le rsultat numrique et rendent le rsultat e e e e e quantitativement inutilisable. Pour amliorer le rsultat il faut diminuer h (donc k) et augmenter le domaine spatial e e discret, on augmente la mmoire et le temps de calcul ncessaires. e e
1.0 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6
t25
t50
t75 t100
PSfrag replacements
0.8 1.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
:v :u
5.5.2
Equation de la chaleur
On se propose de rsoudre numriquement lquation de la chaleur dans un ouvert de R2 : e e e 2 2 u (x, y, t) = d u (x, y, t) + u (x, y, t) t x2 y 2 u(x, y, 0) = (x, y) u(x, y, t) = 0 pour (x, y, t) ]0, a[2 ]0, T ] pour (x, y) ]0, a[2 pour (x, y, t) (]0, a[2 ) [0, T ] ,
Comme dans le cas des fonctions a une dimension spatiale, on utilise la formule de Taylor ` pour approximer les drives partielles de u . Les valeurs de v sur la grille discr`te sont e e e n notes vm,p , o` n N et m, p Z . Pour simplier on suppose que (xm , yp ) = (mh, ph) , e u cest-`-dire que lon a une grille spatiale carre. a e Le schma de Crank-Nicolson scrit e e
n+1 n vm,p vm,p d = k 2 n n n n n n vm+1,p 2vm,p + vm1,p vm,p+1 2vm,p + vm,p1 + h2 h2 n+1 n+1 n+1 n+1 n+1 n+1 vm+1,p 2vm,p + vm1,p vm,p+1 2vm,p + vm,p1 + + h2 h2
5.5. APPLICATIONS
89
Souvent on utilise les masques pour exprimer la relation entre v n et v n+1 : 0 d 0 0 d 0 d 2(1 + 2d) d v n+1 = d 2(1 2d) d v n . 0 d 0 0 d 0
(5.10)
Cette relation est valable pour tous les points intrieurs au domaine spatial, i.e. les points e 2 (xm , yp ) ]0, a[ . Les valeurs sur le bord tant imposes par le probl`me, il faut dterminer e e e e n+1 vm,p en tout point intrieur, grce a la relation ci-dessus. Pour cela il faut rsoudre un e a ` e syst`me linaire que lon va construire maintenant. e e Supposons que [0, a] est discrtis en M + 2 points, les points intrieurs sont alors dans e e e lensemble {(xm , yp ) / 1 m, p M } , car x0 = y0 = 0 et xM +1 = yM +1 = a . En utilisant lordre lexicographique, i = m+(p1)M , pour ordonner les points intrieurs, e on peut crire la relation (5.10) sous forme matricielle : e A V n+1 = B n , o` B n RM et, pour i = 1, . . . , M 2 : u
n n n n n Bin = 2(1 2d) vm,p + d (vm+1,p + vm1,p + vm,p+1 + vm,p1 ). n+1 On a vm,p = Vin+1 et, pour les quatre voisins, n+1 n+1 vm1,p = Vi1 , n+1 n+1 vm+1,p = Vi+1 , n+1 n+1 vm,p1 = ViM , n+1 n+1 vm,p+1 = Vi+M .
2
(5.11)
Tous les lments de la ie ligne, 1 i M 2 , de la matrice A sont nuls, sauf Ai,i = et, ee si si i MN + 1 : i MN 1 : Ai,i1 = ; si Ai,i+1 = ; si i>M : i < M2 M : Ai,iM = ; Ai,i+M = ;
90
La matrice A est tridiagonale par blocs, chacun des M 2 blocs est carr dordre M et A e 2 2 est une matrice (M M ) . Pour chaque pas de temps tn il faut rsoudre (5.11), cest-`-dire un syst`me dordre M 2 . e a e On montre que la matrice A est dnie positive et admet une dcomposition de Cholesky : e e t A = L L , o` L est une matrice triangulaire infrieure. Cette dcomposition simplie la u e e rsolution du syst`me (5.11). An dconomiser de la place mmoire, on utilise de plus e e e e une reprsentation adapte a la matrice A qui est une matrice creuse (angl. sparse). e e ` On propose le code Scilab suivant :
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 // Param`tres du probl`me e e d=0.5; // coefficient de diffusion a=%pi; // domaine spatial ]0,a[x]0,a[ T=1; // domaine temporel [0,T] // Augmenter la mmoire de Scilab e stacksize(5000000); // Param`tres du schma e e h=0.05; MM=ceil(a./h)+1; M=MM-2; xx=linspace(0,a,MM); x=xx(1,2:MM-1); k=0.01; N=ceil(T/k)+1; t=linspace(0,T,N); mu=k./h^2; M2=M^2; v=zeros(MM,MM,N);
// // // // //
pas de discrtisation spatiale e nb. de points dans [0,a] nb. de points dans ]0,a[ intervalle discret despace de 0 a a ` intervalle discret despace sans 0 et a
mprintf(d=%f\ta=%f\tT=%f\n,d,a,T) mprintf(h=%f\tk=%f\n,h,k) mprintf(mu=%f\tlambda=%f\n,mu,k/h) mprintf(N=%i\tM=%i\tM^2=%i\n,N,M,M2) // Condition initiale a evaluer sur les points intrieurs ` e deff([z]=phi(x,y),z=sin(x).*sin(y)); v(2:MM-1,2:MM-1,1)=eval3d(phi,x,x); // Construction des matrices A et B // noter que lon nutilise que des matrices creuses D1=sparse([1:M2-1;2:M2],ones(M2-1,1),[M2 M2]).. -sparse([M:M:M2-M;1+M:M:M2-M+1],ones(M-1,1),[M2 M2]); DM=sparse([1:M2-M;1+M:M2],ones(M2-M,1),[M2 M2]); A=-(d*mu).*(D1+DM); // triangle sup. de A A=A+A; // A est symtrique e B=-A; A=A+2.*(1+2*d*mu).*speye(M2,M2); // on ajoute la diagonale B=B+2.*(1-2*d*mu).*speye(M2,M2);
5.5. APPLICATIONS
91
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
clear D1,DM; mprintf(Stacksize=%i, used=%i\n,stacksize()) // Factorisation de A par la mthode de CHOLESKY e spchoA=chfact(A); clear A; // Affichage de la condition initiale xbasc(); rect=[0,a,0,a,0,1]; xtitle([Condition initiale]); plot3d(xx,xx,v(:,:,1),233,81,"x@y@v(x,y,0)",[2,1,4],rect); // Boucle sur le temps t_n : a chaque fois on rsout le syst`me ` e e // A V^(n+1) = B V^(n) gr^ce a la dcomposition de CHOLESKY a ` e // on obtient le vecteur V=A^{-1} B V^(n) // ensuite reconstruction de v(x,y,t(i)) et affichage for i=2:N V=chsolve(spchoA,B*matrix(v(2:MM-1,2:MM-1,i-1),M^2,1)); v(2:MM-1,2:MM-1,i)=matrix(V,M,M); xclear(); xtitle([Solution a linstant t(+string(i)+)=+string(t(i))]); ` plot3d(xx,xx,v(:,:,i),233,81,"x@y@u(x,y,t)",[2,1,4],rect) end clear B,V;
On obtient lachage des param`tres : e d=0.500000 a=3.141593 T=1.000000 h=0.050000 k=0.010000 mu=4.000000 lambda=0.200000 N=101 M=62 M^2=3844 Stacksize=5000000, used=1024430 et, les graphes suivants :
Condition initiale (x, y) = sin(x) sin(y)
u(x, y, 0)
1.0
0.5
PSfrag replacements
0.0
u(x, y, 1) /2 y
0 0
/2 x
92
u(x, y, 1)
1.0
0.5
PSfrag replacements
0.0
u(x, y, 0) /2 y
0 0
/2 x
On peut de nouveau facilement calculer la solution u de lquation aux drives partielles e e e et comparer a chaque tape la solution discr`te v a u. On utilise pour cel` ` e e ` a err (tn ) = sup |u(xm , yp , tn ) v(xm , yp , tn )|
(xm ,yp )
et
err2 (tn ) = h
(xm ,yp )
1/2
err2
2 103
err
PSfrag replacements
4 104
0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Pour voir leet de la diusion sur une donne initiale irrguli`re, on se propose dtudier e e e e 2 sur ]0, 1[ lvolution de la fonction (x, y) = I[1/3,2/3] (x, y) . On obtient des rsultats nue e mriques satisfaisants pour = k/h petit. e En adaptant le code prcdent, on obtient les rsultats suivants : e e e d=0.010000 a=1.000000 T=1.000000 h=0.030000 k=0.002000 mu=2.222222 lambda=0.066667 N=501 M=33 M^2=1089 Stacksize=5000000, used=672021
5.5. APPLICATIONS
93
Solution a linstant t = 0.1 `
u(x, y, 0)
1
u(x, y, 0.1)
1
0.33
1
1 1.0 0.5 0.5
1.0
0.67 1
1.0
PSfrag replacements
PSfrag replacements
1.0 0.5
0.5
Solution a linstant t = 1 `
u(x, y, 0.5)
1
u(x, y, 1)
1
1 1.0
PSfrag replacements
PSfrag replacements
y
0 0
0.5
Sur ces graphiques on a reprsent aussi quelques lignes de niveau de v , ceci an de e e c e visualiser la diusion du cube initial. Une autre faon de reprsenter v est dassocier des niveaux de gris aux valeurs que prend la fonction, ce qui donne les gures suivantes, o` 0=blanc et 1=noir : u
v(xm , yp , 0)
v(xm , yp , 1)
Remarque : Pour les images des pages 9 et suivantes, on a 255=blanc et 0=noir ; on y a utilis des schmas moins coteux. En eet, lalgorithme de cette section qui utilise la e e u mthode de Crank-Nicolson, est tr`s prcis, mais il ncessite beaucoup de place mmoire e e e e e et de temps de calcul. En dimension spatiale deux on utilise souvent dautres mthodes. e
94
5.5.3
Syst`me de raction-diusion e e
Un exemple de syst`me de raction-diusion qui engendre des rayures a t propos par e e ee e Meinhardt 2 : 2 c1 = d c + c1 s2 c c 1 c 1 t r c 2 c2 s1 = dc c2 + 2 c c2 t r s1 t = ds s1 + s (c1 s1 ) s 2 t = ds s2 + s (c2 s2 ) r = c 2 s2 + c 2 s1 r r 1 2 t e e o` les cinq morphog`nes c1 , c2 , s1 , s2 et r sont dnies sur [0, M ]2 R+ . Le comu portement du syst`me dpend des param`tres rels dc , ds , c , s et r et des conditions e e e e initiales. Pour ce type de simulation on a besoin dun nombre lev ditrations, N , et il devient e e e ncessaire dutiliser un programme compil. Ainsi le syst`me ci-dessus a t programm e e e ee e en langage C en utilisant les librairies de Megawave3 . On a utilis un schma explicite, e e progressif en temps et centr en espace, avec des conditions au bord priodiques. Les e e conditions initiales sont donnes par e r r (1 + 1 (x, y)) , c2 (x, y, 0) = (1 + 2 (x, y)) , c1 (x, y, 0) = 2c 2c r 2 s1 (x, y, 0) = s2 (x, y, 0) = , r(x, y, 0) = r3 , 2c 4c pour (x, y) [0, M ]2 et o` 1 et 2 est un bruit uniforme, de valeurs dans [ , ] , avec u > 0 donn. La discrtisation spatiale est xe a h = 1 dans les deux directions, le pas e e e ` de temps k est obtenu en choisissant : k = h2 = . Le choix des param`tres pour les e premi`res simulations est indiqu dans la table suivante : e e = 0, 0100 M = 100 N = 30.000 dc = 0, 0500 ds = 0, 1000 = 0, 1000 c = 0, 4000 s = 0, 4000 r = 0, 6000
Models of Biological Pattern Formation, Academic Press, 1982 c logiciel disponible a http ://www.cmla.ens-cachan.fr/Cmla/Megawave `
5.5. APPLICATIONS
95
96
Pour les param`tres du tableau suivant on obtient les zbrures reprsentes en dessous. e e e e = 0, 0050 M = 100 N = 50.000 dc = 0, 1500 ds = 0, 2000 = 0, 1000 c = 0, 4000 s = 0, 4000 r = 0, 8000
c1 (xm , yp, kN )
c2 (xm , yp , kN )
s1 (xm , yp , kN )
s2 (xm , yp , kN )
r(xm , yp , kN )
5.5. APPLICATIONS
Annexes
x1 . Un point x = (x1 , . . . , xd ) de Rd considr comme vecteur sera not x = . ee e . xd dont la transpos , xt = ( x1 xd ) , est une matrice 1 d . e
d 1/2
x2 i
i=1
x et y scrit x.y = xt y = e
i=1
xi y i .
f (a) = fxi (a) xi la drive partielle de f en a par rapport a la variable xi , i = 1, . . . , d . e e ` k f Pour k N on a Dk f (a) = (a) . i xk i Soit f C1 (Rd , R) , pour a = (a1 , . . . , ad ) Rd , on note Di f (a) = On note Df (a) = D1 f (a) D2 f (a) . . . Dd f (a) cation direntielle dfa L(Rd , R) . e Le vecteur gradient est le vecteur colonne (d 1) la matrice ligne (1 d) associ a lapplie` f (a) = Df (a)t .
Dans la suite on va souvent utiliser la notation des multi-indices de L. Schwartz. Soit = (1 , . . . , d ) Nd et x = (x1 , . . . , xd ) Rd , on pose || = 1 + + d , Pour f C|| (Rd , R) , on note D f (a) = D1 D2 . . . Dd f (a) = 1 2 d || f (a) , x1 xd 1 d ! = 1 ! 2 ! d ! et x = x1 x2 xd . 2 1 d
100
Soit F C1 (Rd , Rn ) et x = (x1 , . . . , xd ), on note D1 F1 (x) D2 F1 (x) . . . D1 F2 (x) D2 F2 (x) . . . DF (x) = . . . . . . D1 Fn (x) D2 Fn (x) . . .
Dd Fn (x)
Cest la matrice (n d) associ a lapplication direntielle dFx L(Rd , Rn ) qui est e ` e compose des d vecteurs colonnes Di F (x). e
A.2
Calcul direntiel e
Soit f C1 (Rd , R), on dnit la drive directionnelle de f au point a Rd et dans la e e e direction de v Rd , par df 1 (a) = lim (f (a + hv) f (a)) . h0 h dv On montre que df (a) = dv
Rd+1
f (a).v
f (a)
PSfrag replacements
Rd
fa a
Soit = (1 , . . . , d ) Nd , la notation des multi-indices de L. Schwartz permet dnone cer facilement les rsultats suivants : e Si P (x) est un polynme de d variables, x1 , . . . , xd , et de degr m, o e alors P (x) =
||m
1 (D P (0)) x . ! m! x , !
(x1 + + xd )m =
||=m
101
Soit un ouvert de Rd et f Cn (, R) . On suppose que 0 et que le segment [0, x] , alors la formule de Taylor avec reste intgral, pour f en 0 R d , scrit e e
n1
f (x) =
p=0 ||=p
1 D f (0) x + !
||=n
n !
1 0
(1 s)n1 D f (sx) ds
x .
! (D f ) (D g) . !!
m! (D f (x + tv)) v . !
A.3
Intgration dans Rd e
On va noter dans cette section dd ou dx la mesure de Lebesgue sur Rd . On rappelle le rsultat suivant. e Thor`me A.3.1 (Changement de variables dans Rd ) e e Soit U et V des ouverts de Rd , un diomorphisme de U sur V et soit f une fonction e intgrable sur V , alors e f (x)dx =
V U
f ((y))|dt(D(y))|dy e
ou encore
U
f dd =
f |dt(D1 )| dd . e
On note Bd (x0 , r) = {x Rd /|x x0 | < r} , la boule ouverte de centre x0 et de rayon r , la boule unit de Rd est note Bd = Bd (0, 1) . e e d1 De mme S (x0 , r) = {x Rd /|x x0 | = r} = Bd (x0 , r) e et S d1 = Bd (0, 1) est la sph`re unit de Rd . e e
x Soit x Rd \ {0} et h(x) = |x| , on obtient un homomorphisme de Rd sur ]0, +[S d1 e 1 en posant (x) = (|x|, h(x)) = (r, s) et (r, s) = rs .
f (x)dd (x) =
Rd ]0,+[S d1
dr
S d1 (0,r)
f (x)dS(x)
Limage de d par est la mesure produit r d1 dr dd1 , o` d1 est une mesure sur u d1 B(S ) .
ANNEXE A. CALCUL DIFFERENTIEL ET INTEGRATION DANS RD
102
Corollaire A.3.1 Soit : R+ R borelienne, on pose f (x) = (|x|). Alors f est intgrable sur Rd si et e d1 seulement si (r)r est intgrable sur R+ . e Corollaire A.3.2 1 Soit R, d N et f (x) = , alors |x| - pour tout < d : f est intgrable sur Bd et nest pas intgrable sur Rd \ Bd ; e e d - pour tout > d : f est intgrable sur R \ Bd et nest pas intgrable sur Bd ; e e - pour = d : f nest intgrable ni sur Bd , ni sur Rd \ Bd . e Corollaire A.3.3 On note d le volume de la boule unit Bd de Rd et d laire de la sph`re unit S d1 de e e e Rd , on a les relations : 2 d/2 d/2 = ( d + 1) d ( d ) 2 2 d/2 ( d ) 2
= d = = rd d
Exemples : On retrouve ainsi, pour d = 2, que le volume dun disque est r 2 (=surface) et que l aire dun cercle est 2r (=prim`tre). e e 4 Pour d = 3 on retrouve 3 = 3 et 3 = 4 .
Proposition A.3.2 Soit f une fonction continue sur Rd , alors pour tout x Rd : f (x) = lim
r0 r d
1 d
f (y) dy = lim
Bd (x,r)
1 d
S d1 (x,r)
r0 r d1
f (y) dSy .
Avec des hypoth`ses plus faibles on a le rsultat suivant : e e Thor`me A.3.2 (Thor`me de drivation de Lebesgue) e e e e e Pour tout f de L1 (Rd ) on a : 1 (a) pour presque tout x Rd : f (x) = lim d f (y) dy ; r0 r d Bd (x,r) 1 (b) pour presque tout x Rd : lim d |f (y) f (x)| dy = 0 . r0 r d Bd (x,r) On dit que x est un point de Lebesgue de f sil vrie (b). e
Dans cette section on va gnraliser la notion de courbe paramtrique dans R d , dnie e e e e d+1 par (X1 (t), . . . , Xd (t)) , et celle de surface dans R , dnie par le graphe dune fonction e z = f (x1 , . . . , xd ) . Dfinition B.1.1 e Soit M Rd et a M , on dit que M est une varit rguli`re de dimension k au voisinage ee e e de a sil existe : U un ouvert de Rk , k d , W un voisinage ouvert de a dans Rd et C1 (U, Rd ) tel que : est un homomorphisme de U sur W M et pour tout u U , est de rang k , e cest-`-dire rang(D(u)) = k pour tout u U . a On dit que est une paramtrisation de M au voisinage de a. e M est une varit rguli`re de dimension k, si M est au voisinage de tout point a une ee e e k-varit rguli`re. ee e e On obtient ainsi une collection de paramtrisation locales et on appelle (W M, ) une e carte et {(W M, )} un atlas de M . En anglais une varit sappelle manifold. ee
Rd W PSfrag replacements a M W
b R
k
Exemples : Soit M une k-varit rguli`re de Rd . Pour k = 1, M est une courbe ; pour ee e e k = 2, cest est une surface ; pour k = d 1, on dit que M est une hypersurface .
104
Thor`me B.1.1 e e Soit M Rd et a = (a1 , . . . , ad ) M , on a quivalence des armations suivantes : e (P) M est une k-varit rguli`re au voisinage de a. ee e e (Z) Il existe W un voisinage ouvert de a dans Rd et C1 (W, Rdk ) tel que : est de rang d k en a, (a) = 0 et 1 (0) = W M . (G) Il existe W un voisinage ouvert de a dans Rd , V un voisinage ouvert de (a1 , . . . , ak ) dans Rk et C1 (V, Rdk ) tel que, si on note = (k+1 , . . . , d ), on a pour tout x W M et pour k + 1 i d : i (x1 , . . . , xk ) = xi . Remarques : a) Avec (P ) on reprsente M localement comme limage dune application : Rk Rd ; e (Z) exprime que M est localement lensemble des points qui annulent une application : Rd Rdk ; nalement (G) montre que M est localement le graphe dune application : Rk Rdk . b) (G) est vrie a une permutation des coordonnes pr`s. e e ` e e Par ailleurs (G) est un cas particulier des reprsentations (P ) et (Z), en eet, si on e se donne , alors (x1 , . . . , xd ) = (x1 , . . . , xk , k+1 (x1 , . . . , xk ), . . . , d (x1 , . . . , xk )) est une paramtrisation, i.e. vrie (P ), et (x1 , . . . , xd ) = (x1 , . . . , xk ) (xk+1 . . . , xd ) e e Rdk est une application qui sannule au voisinage de a , c.`.d. vriant (Z) . a e c) Dans (Z), lhypoth`se que est de rang d k en a, i.e. rang(D(a)) = d k, est e quivalente a demander que da est surjective de Rd sur Rdk , cest-`-dire que a nest e ` a pas un point critique de . d) En utilisant (G) on peut toujours aplatir la k-varit rguli`re M au voisinage de a : ee e e d d Soit : R R deni par (x) = 1 (x), . . . , d (x) o`, pour x = (x1 , . . . , xd ) Rd , u D k+1 et D = . . . Dd Ik 0k,dk Idk .
i (x) =
xi si 1 i k xi i (x1 , . . . , xk ) si k + 1 i d
Rk (a1 , . . . , ak , 0, . . . , 0)
(M )
105
Exemples : 1. Les varits de dimension k = 0 sont les points discrets de Rd . ee Les varits de dimension k = d sont les ouverts de Rd ee Tout sous-espace linaire ou ane de Rd est une varit rguli`re. e ee e e 2. Lensemble {(x1 , x2 ) R2 / x1 x2 = 0} nest pas une varit de R2 : ee au voisinage de (0, 0) lensemble nest pas homomorphe a un segment de droite. e ` 3. Soit M une 2-varit, i.e. une surface, de R3 au voisinage de a = (a1 , a2 , a3 ) . ee Dapr`s (P ) il existe une paramtrisation : U R2 R3 , telle que rang(D) = 2 e e en tout point de U . Si lon note (u, v) = (1 (u, v), 2 (u, v), 3 (u, v)) pour (u, v) R2 on a 1 1 u v 2 2 D = u v 3 3 u v qui est de rang deux si et seulement si lun au moins de ses trois sous-dterminants e est non nul : 1 (1 , 2 ) = u 2 (u, v) u 1 2 v , (2 , 3 ) = u 2 3 (u, v) v u 2 3 v , (3 , 1 ) = u 3 1 (u, v) v u 3 v . 1 v
Dapr`s (G) il existe tel que x3 = (x1 , x2 ) pour tout x = (x1 , x2 , x3 ) de M e Bd (a, ) . On en dduit une paramtrisation en posant 1 (u, v) = u, 2 (u, v) = v et 3 (u, v) = e e (u, v) pour (u, v) dans un voisinage de (a1 , a2 ) . Si lon pose (x1 , x2 , x3 ) = (x1 , x2 ) x3 on a (x1 , x2 , x3 ) 1 (0) pour tout x M Bd (a, ) . Application : Pour S 2 on obtient une paramtrisation au voisinage du ple nord e o a = (0, 0, 1), avec (x1 , x2 ) = 1 x2 x2 on a : 1 2 (u, v) = (u, v, (u, v)) et (x1 , x2 , x3 ) = (x1 , x2 ) x3 . Que devient S 2 si on laplatit au voisinage du ple nord ? o
106
e recouvrements ventuels des cartes sont direntiables , au sens de la e Proposition B.1.1 Si on a deux paramtrisations C1 , et , au voisinage de a M , alors il existe T un e diomorphisme de Rk sur Rk qui permet de passer de lune a lautre. On dnit : e ` e T = 1 : 1 (U ) (U ) 1 (U ) (U ) , pour tout , tel que (U ) (U ) = .
Rd a
PSfrag replacements
U U Rk T
B.2
Dfinition B.2.1 e On appelle espace tangent a la k-varit rguli`re M au point a M et on note T a M , ` ee e e le sous-espace vectoriel de dimension k, image de Rk par lapplication linaire db , o` e u a = (b) . Soit a M , a = (b), alors si x a + Ta M , il existe v Rk tel que ax = x a = db (v). 1 1 u1 uk . . On a Db = . , on en dduit immdiatement : e e . . . d d u1 uk |u=b
k
Equation paramtrique de a + Ta M : xi = ai + e
j=1
i (b) vj pour i = 1, . . . , d . uj
Cas dune hypersurface : 1 1 u1 ud1 . . . = 0. Si k = d1 on obtient lquation cartsienne de a+Ta M : e e . . . . . . d d (xd ad ) u1 ud1 (x1 a1 )
107
U Rk
Comme M est une hypersurface de Rd , on peut reprsenter M au voisinage de a grce e a au graphe de (u1 , . . . , ud1 ) , avec a = b, (b) et b = (a1 , . . . , ad1 ). On a Db = = ( t1 0 0 1 u1 u2 ud1 1 0 . . . 0 1 . .. . . . 0 0 . . .
td1 )
et la famille {t1 , . . . , td1 } est une base de lhyperplan Ta M . On en dduit le vecteur normal unit a M en a : n(a) = e e`
u1 . . 1 . . 1 + | |2 ud1 1
Exemple : On suppose que M est une surface dans R3 , reprsente localement par le e e graphe dune fonction . Avec les notations prcdentes on a : a = (a1 , a2 , a3 ) M , e e (a1 , a2 ) = a3 et b = (a1 , a2 ), et une paramtrisation locale est donne par (u, v) = e e (u, v, (u, v)) . 1 0 1 et lquation cartsienne du plan tangent ane scrit : e e e Donc Db = 0 (b) (b) u v (x1 a1 ) (a1 , a2 ) + (x2 a2 ) (a1 , a2 ) (x3 a3 ) = 0 . u v
n(a) =
2 v
u v
108
u . pour u2 +v 2 < 1 et (u, v) = 1 u2 v 2 , on obtient n(u, v, (u, v)) = v 1 u2 v 2 On va sintresser maintenant a lorientation de M . Rappelons dabord quun rep`re direct e ` e de Rk vrie dt(u1 , . . . , uk ) > 0 et que lon dnit ainsi une orientation sur Rk . e e e Soit une paramtrisation locale de la k-varit rguli`re M au voisinage de a = (b) e ee e e M . Limage par d (b) du rep`re direct est une base de Ta M . On dnit ainsi un rep`re e e e donc une orientation sur Ta M . Soit une autre paramtrisation locale de M au voisinage de a. On a le changement de e 1 param`tres T (u) = ( )(u) pour u 1 ( (U ) (U )) et d = d dT . e On obtient la mme orientation pour Ta M si dt(DT (u)) > 0 . e e Dfinition B.2.2 e e On dit que la varit rguli`re (M, { } ) est orient si pour tout = : dt(DT ) > 0 ee e e e en tout point u o` T est dni. u e Une varit M est orientable sil existe des paramtrisations locales qui orientent M . ee e Remarque : Une varit nest pas toujours orientable, la bande de Mbius est un contreee o exemple. Exemples : Courbes dans Rd Soit : R Rd une paramtrisation dune courbe (k = 1), est le vecteur tangent e a la courbe, on lutilise pour dnir une orientation sur la courbe. ` e Hypersurfaces dans Rd Plus haut on a calcul le vecteur normal unit n(a) a une hypersurface de Rd , il dnit e e ` e une orientation de M . En particulier si M = , o` est un ouvert born de Rd , on pourra toujours choisir u e n(a) pointant vers lextrieur de . e
109
PSfrag replacements
Bande de Mbius (1858) o Soit R > 1, on consid`re une ne tige de longueur 2 dans le plan xz : e (x, z) = (R + t, 0) pour 1 < t < 1 . On fait tourner la tige de 180 degrs autour de son centre et, en mme temps, on lui e e fait faire un tour complet autour de laxe Oz , la surface ainsi obtenue est une bande de Mbius. o Montrer quune paramtrisation de la bande de Mbius est donne par e o e (R + t cos ) cos 2 (t, ) = (R + t cos ) sin 2 pour 1 < t < 1 et 0 . t sin Vrier que la bande de Mbius nest pas orientable. e o
z
PSfrag replacements x y
110
Hyperbolo ` une nappe. de a Soit M = (x, y, z) R3 / x2 + y 2 z 2 = 1 . Montrer que M est une 2-varit rguli`re de R3 . Donner une paramtrisation de M . ee e e e
z
PSfrag replacements
Cne. o Soit M = (x, y, z) R3 / x2 + y 2 z 2 = 0 . Donner une paramtrisation de M . Est-ce que M est une 2-varit rguli`re de R 3 ? e ee e e
z
x PSfrag replacements
111
B.3
Intgrale de surface e
Dans cette section on va introduire lintgrale dune fonction dnie sur une varit et e e ee mesurer laire de portions de varits. On parle dintgrale de surface et aire de surface ee e mme si k = 2 . e Avant de passer a laire des varits de Rd on va sintresser au calcul, plus simple, du ` ee e d volume dun paralllpip`de de R . ee e Proposition B.3.1 Soit P le paralllpip`de engendr par les k vecteurs l.i. a1 , . . . , ak de Rd : ee e e
k
P = On note A = ( a1
xR
x=
i=1
ti ai , ti [0, 1]
. . . ak ) la matrice d k de vecteurs colonnes a1 , . . . , ak , alors si k = d : vol(P ) = dt(At A)1/2 , e si k < d : aire(P ) = dt(At A)1/2 = volk (P ) . e
t a1 a1i . t . Remarque : On a A = ( a1 . . . ak ) avec ai = . et A = . , . . adi at k a1 .a1 a1 .a2 . . . a1 .ak a2 .a1 a2 .a2 . . . a2 .ak donc At A = . e e . . est une matrice carre symtrique k k .. . . . . . . . ak .a1 ak .a2 . . . ak .ak
a1i . 2. Dans R3 on consid`re deux vecteurs l.i. ai = . (i = 1, 2) qui engendrent le e . adi paralllpip`de P . ee e Montrer que aire(P ) = |a1 a2 | .
1. Dans 3 on consid`re le paralllpip` engendr par : R e e e ede e a11 0 a13 a1 = 0 , a2 = a22 et a3 = a23 , reprsenter P . e 0 0 a33 Montrer que vol(P ) = |a11 a22 a33 | = base hauteur .
112
Soit S une varit rguli`re compacte de dimension k, cest-`-dire tel quil existe une kee e e a varit rguli`re M avec S M et S compact. Soit : Q Rk S une paramtrisation ee e e e de S que lon suppose tre une bijection de Int(Q) sur S . e
l
Qi ,
k
PSfrag replacements
Qi
Rk
aire((Qi ))
On va procder par approximations : localement on remplace S par son plan tangent de e faon a avoir aire((Qi )) aire(dqi (Qi )) . c ` Si on note {e1 , . . . , ek } la base canonique de Rk , le cube Qi est engendr par les vecteurs e {hej }j=1,...,k et dqi (Qi ) est le paralllpip`de engendr par les vecteurs h ee e e (qi ) . uj j=1,...,k
Posons Ai = ( h u1 (qi ) h u2 (qi ) h uk (qi ) ) = h D(qi ) , on a
1/2
Donc aire(S)
dt Dt i Dqi e q
i=1
1/2
volk (Qi ) .
dt Dt Du e u
1/2
du
f ((u)) dt Dt Du e u
1/2
du
113
Soient et deux paramtrisations de S : e k Rk S et (Q) = (Q) . :QR S ,:Q Il existe un diomorphisme T de Q sur Q : T = 1 et dt(DT ) = 0 sur Q . e e Si lon note u = T(u) et comme du = du dTu on a e dt Dt Du e u donc dt Dt Du e u
1/2 1/2 t = |dt(DTu )| dt Du Du e e e e
1/2
du =
Q
dt Dt (u) DT (u) e T
1/2
|dt(DTu )| du = e
e Q
t dt Du Du e e e
1/2
du .
Do` le rsultat suivant : u e Proposition B.3.2 La dnition de laire dune varit rguli`re compacte, respectivement de lintgrale sur e ee e e e une telle varit, est indpendante de la paramtrisation choisie. ee e e
Exemples : 1. Soit C une courbe rguli`re (k = 1) dans Rd , et : [a, b] R Rd une parame e e trisation. Comme Du = (u) , on a Dt Du = | (u)|2 et, en remarquant quici laire est u en fait la longueur :
b b d
long(C) =
a
| (u)| du =
i (u)2 du
i=1
2. Soit S une surface rguli`re (k = 2) dans Rd , et une paramtrisation e e e : Q R2 S . (u, v) (u, v) = (1 (u, v), . . . , d (u, v)) On a D = u v
d
do` Dt D = u
2
. u u . u v d
. u v . v v
et, si on note
d
E= . = u u on obtient
i=1
i u
,F = . = u v
i=1
i i et G = . = u v v v
i=1
i v
aire(S) =
Q
EG F 2 dudv
114
3. Soit S une hypersurface rguli`re (k = d 1) dans Rd , et une paramtrisation e e e : Q Rd1 S , (u1 , . . . , ud1 ) (u1 , . . . , ud1 , (u1 , . . . , ud1 ))
0 1 . . . 0
u2
.. .
t
0 0 . . . 1
ud1
et on montre que
d1 i=1
dt D D = e do` u aire(S) =
Q
ui
= 1 + | |2
1 + | (u)|2 du
Cas particulier : d = 3 et k = 2, si S est reprsente comme graphe de : z = (u, v) e e pour (u, v) Q R2 on trouve aire(S) =
Q
1+
dudv .
U (x).n(x) dSx
PSfrag replacements
u (x) dx = xi
o` ni est la i-`me coordonne du vecteur normal unit extrieur. u e e e e Soient u et v deux champs scalaires de classe C1 sur , alors, pour tout 1 i d , on a la formule dintgration par parties dans Rd e u (x)v(x) dx = xi u(x)v(x)ni (x) dSx u(x)
v (x) dx xi
116
Exercice : crire cette formule dans le cas d = 1 , avec = [a, b] . e Si on suppose que u et v sont de classe C2 sur , on dduit de ce qui prc`de les e e e formules de Green (1828) :
u(x) v(x) dx =
u(x). v(x) dx
(I)
u(x) v(x) dx =
u(x) v(x) dx +
Remarque : Les rsultats prcdents sobtiennent a partir de la formule de Stokes qui e e e ` est un rsultat de ltude des formes direntielles. En particulier on montre que si e e e est une varit orient et rguli`re par morceaux les noncs prcdents restent vrais. ee e e e e e e e
D(x0 , )
u(x, t) dx .
La variation de masse m (t) , ngative si plus de particules sortent de quil nen rentrent, e est value de deux faons direntes. e e c e
117
Mthode 1. Soit P un paralllpip`de rectangle innitsimal, centr en p = (p1 , p2 , p3 ), e ee e e e de faces parall`les aux axes de coordonnes et de cts x1 , x2 et x3 . On note e e oe x xi 3 Fi = x R / xi = pi 2 et pour j = i : |xj pj | 2 j les six faces de P . Pour obtenir la variation de masse on va estimer le nombre de particules qui sortent du paralllpip`de P . ee e
PSfrag replacements x3 x2
F1
x1
x3 p
+ F1
x1
x2
+ On peut supposer que le ux de particules au centre de la face F1 est approximativement + 1 U1 donne par U1 (p) + 2 x1 (p)x1 et que la quantit de particules traversant F1 par unit e e e
1 U1 (p)x1 2 x1
x2 x3 x2 x3
(1) .
Finalement la variation du nombre de particules dans P , due a des dplacements dans la ` e U1 direction x1 , est (2) (1) = x1 (p)x1 x2 x3 . Cette valeur est positive si plus de particules rentrent dans P que nen sortent par les + faces F1 et F1 .
On obtient des rsultats analogues pour la direction x2 (faces F2 ) et x3 (faces F3 ). e La variation par unit de temps du nombre de particules dans P est donc donne par e e U2 U3 U1 x1 (p) + x2 (p) + x3 (p) x1 x2 x3 = div(U (p))x1 x2 x3 .
U1 (p) x1 x1 2 U1 x1 x1 (p) 2
En faisant la somme sur tous les paralllpip`des rectangles innitsimaux P contenus ee e e dans on a m (t) = div(U (x)) dx .
Mthode 2. On va mesurer maintenant le nombre de particules qui traversent au cours e du temps. Soit a et Sa un lment de surface susamment petit pour tre ee e proche de a + Ta () et soit n(a) le vecteur normal unit extrieur a en a . e e ` Alors la variation par unit de temps du nombre de particules a travers Sa est donne e ` e par U (a, t).n(a) aire(Sa ) = u(a, t) aire(Sa ) V (a, t).n(a) , et nalement en faisant la somme sur tous les lments de surface Sa on obtient comme valuation de la variation de ee e masse m (t) = U (x, t).n(x) dSx .
Conclusion : en identiant les variations m (t) obtenus par les deux mthodes on trouve e lnonc du thor`me de la divergence. e e e e
118
Le thor`me de la divergence permet de donner linterprtation suivante de la divergence e e e dun champ de vecteurs U (x) : div(U (x0 )) = lim 1 d
d S d1 (x0 , )
U (x).n(x) dSx ,
cest-`-dire que la divergence de U en x0 est le taux de particules (par unit de volume) a e qui divergent de x0 .
u t
dx = 0 .
Si on a lgalit ci-dessus pour tout volume on en dduit la loi dquilibre du milieu e e e e u u + div(U ) = + t t u.V + u div(V ) = 0 .
On dit que le ux est incompressible si div(U ) = 0 , ici u = C te et div(V ) = 0 . Exemple : On consid`re un tuyau T qui se resserre (cf. gure) et dans lequel coule un e liquide, si le liquide est incompressible on a
div(U ) dx =
U.n dS = 0 or est
U.n dS =
U.n dS ,
S2
S1
S2
119
On dnit le vecteur rotationnel de U en x = (x1 , x2 , x3 ) par : e U3 U2 ( x2 x3 )(x) 3 1 rot(U (x)) = ( U3 U1 )(x) . x x 1 2 ( U1 U2 )(x) x x Soit S une surface (2-varit) orient et compacte dans R3 tel que S est une courbe ee e (1-varit) orient et rguli`re. ee e e e
n(a) t(a) S t(a) n(a)
PSfrag replacements M
Quitte a inclure S dans une 2-varit orient rguli`re M , on dnit en tout point x de ` ee e e e e S M le vecteur unit normal a S , n(x) , pour a S M on dnit en plus le vecteur e ` e unit tangent a S , t(a) . e ` On choisit t(a) tel que S se trouve a gauche en parcourant S, cest-`-dire que t(a) n(a) ` a est le vecteur unit normal a S et extrieur a S (cf. gure), alors, grce au thor`me de e ` e ` a e e
Stokes
U (x).t(x) dx
Soit n S 2 et D(x0 , ) le disque de centre x0 et de rayon , perpendiculaire a n , alors ` U (x).t(x) dx = rot(U (x)).n dSx et donc 1 2
D(x0 , ) D(x0 , )
U (x).t(x) dx .
D(x0 , )
La composante, dans la direction n , du rotationnel de U , en x0 , est le taux de particules (par unit de surface) qui tournent autour de x0 par rapport a laxe de rotation n . e `
PSfrag replacements rot(U (x0 )) n x0 D(x0 , )
Si rot(U ) = 0 on dit que U est un champ irrotationnel. Remarque : Si u est un champ scalaire sur Rd on a u(x0 ) = lim 1 d
d S d1 (x0 , ) 0
u(x)n(x) dSx .
x1 x2 x3
u=u
et le gradient du champ de vecteurs U dans la direction v qui est dni par e v. U1 (x0 ) U (x0 + v) U (x0 ) (v )U (x0 ) = lim = v. U2 (x0 ) = (v. )U (x0 ) = DU (x0 ) v . 0 v. U3 (x0 )
Des notions utiles pour crire des syst`mes dquations aux drives partielles sont e e e e e le Laplacien du champ de vecteurs U qui est dni par e U1 U = div(U ) rot rot(U ) = ( . )U = U2 ; U3