Analyse Numérique des Équations aux Dérivées Partielles PDF - A.U. 2022/2023

Analyse numérique des équations aux dérivées partielles Méthode des différences finies Jawad SALHI Faculté des Sciences et Techniques - Errachidia A.U. 2022/2023 Jawad SAL...

Analyse numérique des équations aux dérivées partielles Méthode des différences finies Jawad SALHI Faculté des Sciences et Techniques - Errachidia A.U. 2022/2023 Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 1 / 82 Plan 1 Introduction Quelques exemples d’équations aux dérivées partielles Notion de problème bien posé Classification des équations aux dérivées partielles 2 Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Analyse de la méthode des différences finies 3 Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Schéma implicite et schéma de Crank-Nicolson 4 Méthode des différences finies pour l’équation des ondes Le schéma explicite naturel Le schéma implicite naturel 5 Notes bibliographiques et remarques Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 2 / 82  Introduction Pour aborder le calcul numérique (à l’aide d’un outil informatique) des solutions d’un problème réel, on passe par les étapes suivantes : Description qualitative des phénomènes physiques : Cette étape, effectuée par des spécialistes des phénomènes que l’on veut quantifier (ingénieurs, chimistes, biologistes etc...) consiste à répertorier tous les mécanismes qui entrent en jeu dans le problème qu’on étudie. Modélisation : Il s’agit, à partir de la description qualitative précédente, d’écrire un modèle mathématique. En d’autres termes, la modélisation mathématique est la science de représenter (ou de transformer) une réalité physique en des modèles abstraits accessibles à l’analyse et au calcul. On supposera ici que ce modèle amène à un système d’équations aux dérivées partielles (EDP). Selon les hypothèses effectuées, la modélisation peut aboutir à plusieurs modèles, plus ou moins complexes. Dans la plupart des cas, on ne saura pas calculer une solution analytique, explicite, du modèle ; on devra faire appel à des techniques de résolution approchée. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 3 / 82  Introduction Analyse mathématique du modèle mathématique : Même si l’on ne sait pas trouver une solution explicite du modèle, il est important d’en étudier les propriétés mathématiques, dans la mesure du possible. Il est bon de se poser les questions suivantes : Le problème est-il bien posé, c’est-à-dire y a-t-il existence et unicité de la solution ? Les propriétés physiques auxquelles on s’attend sont elles satisfaites par les solutions du modèle mathématique ? Si l’inconnue est une concentration, par exemple, peut-on prouver que la solution du modèle censé représenter le modèle physique est toujours positive ? Y a-t-il continuité de la solution par rapport aux données ? Discrétisation et résolution numérique : Un problème posé sur un domaine continu (espace-temps) n’est pas résoluble tel quel par un ordinateur, qui ne peut traiter qu’un nombre fini d’inconnues. Pour se ramener à un problème en dimension finie, on discrétise l’espace et/ou le temps. Si le problème original est linéaire on obtient un système linéaire. Si le problème original est non linéaire on aura un système non linéaire à résoudre. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 4 / 82  Introduction Analyse numérique : Il s’agit maintenant de l’analyse mathématique du schéma numérique. En effet, une fois le problème discret obtenu, il est raisonnable de se demander si la solution de ce problème est proche et en quel sens, du problème continu. De même, si on doit mettre en oeuvre une méthode itérative pour le traitement des non-linéarités, il faut étudier la convergence de la méthode itérative proposée. Mise en oeuvre, programmation et analyse des résultats : La partie mise en oeuvre est une grosse consommatrice de temps. Actuellement, de nombreux codes existent, qui permettent en théorie de résoudre "tous" les problèmes. Il faut cependant procéder à une analyse critique des résultats obtenus par ces codes, qui ne sont pas toujours compatibles avec les propriétés physiques attendues... Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 5 / 82  Introduction Quelques exemples d’équations aux dérivées partielles Équation de Poisson Soit d la dimension de l’espace (en général, d = 1, 2 ou 3), l’équation de Poisson s’écrit : − div(κ∇u) = f , En une dimension d’espace, l’équation de Poisson s’écrit donc : −(κu 0 )0 = f. Le réel κ est un coefficient de diffusion et le terme −κ∇u est un flux de diffusion. Avec un coefficient constant κ = 1, l’équation de Poisson s’écrit en une dimension d’espace −u 00 = f. En deux dimensions d’espace, elle s’écrit −∆u = f avec 2 2 2 2 ∆u = ∂xx u + ∂yy u, où ∂xx u (respectivement ∂yy u) désigne la dérivée partielle seconde de u par rapport à x (respectivement y ). Dans le cas f = 0, on obtient l’équation de Laplace : −∆u = 0. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 6 / 82  Introduction Quelques exemples d’équations aux dérivées partielles Équation de la chaleur On note x la variable d’espace, c’est-à-dire un point de Ω ⊂ Rd (d = 1, 2 ou 3), et t la variable de temps. Dans Ω les sources de chaleur sont représentées par une fonction f (t, x ), tandis que la température est une fonction inconnue u(t, x ). En combinant la loi de conservation de l’énergie et la loi de Fourier, on obtient une équation de la temérature u : ut − ∆u = f. Il faut ajouter à cette équation qui est valable dans tout le domaine Ω, une relation, appelée condition aux limites, qui indique ce qui se passe à la frontière ou au bords ∂Ω du domaine, et une autre relation qui indique quel est l’état initial de la température. Par convention, on choisit l’instant t = 0 pour être le temps initial, et on impose une condition initiale u(t = 0, x ) = u0 , où u0 est la fonction de distribution initiale de température dans le domaine Ω. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 7 / 82  Introduction Quelques exemples d’équations aux dérivées partielles En ce qui concerne la condition aux limites, cela dépend du contexte physique. Si le domaine est supposé baigner dans un thermostat à température nulle, alors la température vérifie la condition aux limites de Dirichlet u(t, x ) = 0 pour tout x ∈ ∂Ω et t > 0. Ainsi, on obtient l’équation de la chaleur :  ut − ∆u = f   pour (x , t) ∈ Ω × R∗+ , u(t, x ) = 0 pour (x , t) ∈ ∂Ω × R∗+ ,  pour x ∈ Ω.  u(t = 0, x ) = u 0 Il s’agit d’une équation d’ordre 1 en temps et d’ordre 2 en espace. Si le domaine est supposé thermiquement isolé de l’extérieur, alors le flux de chaleur sortant au bord est nul et la température vérifie la condition aux limites de Neumann ∂u (t, x ) := n(x ) · ∇u(x ) = 0 pour tout x ∈ ∂Ω et t > 0, ∂n où n est la normale extérieur unité de Ω. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 8 / 82  Introduction Quelques exemples d’équations aux dérivées partielles Équation des ondes L’équation des ondes modélise des phynomènes de propagation d’ondes ou de vibration. Par exemple, en deux dimensions d’espace elle est un modèle pour étudier les vibrations d’une membrane élastique tendu (comme la peau d’un tambour). Au repos, la membrane occupe un domain plan Ω. Sous l’action d’une force normale à ce plan d’intensité f , elle se déforme et son déplacement normal est noté u. On suppose qu’elle est fixée sur son bord, ce qui donne une condition aux limites de Dirichlet. L’équation des ondes dont u est solution est donnée par :    utt − ∆u = f dans Ω × R∗+ , sur ∂Ω × R∗+ ,  u = 0    u(t = 0) = u0 dans Ω,  ut (t = 0) = u1 dans Ω.   Remarquons qu’il s’agit d’une équation du deuxième ordre en temps et qu’il faut donc deux conditions initiales pour u. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 9 / 82  Introduction Notion de problème bien posé Le mathématicien Jacques Hadamard a donné une définition de ce qu’est un "bon" modèle, en parlant de problème bien posé. On décide de noter f les données (le second membre, les données initiales, le domaine etc), u la solution recherchée et A l’opérateur qui agit sur u. Il s’agit ici de notations abstraites, A désignant à la fois l’équation aux dérivées partielles et le type de conditions initiales ou aux limites. Le problème est donc de trouver u solution de A(u) = f. (1) Definition On dit que le problème (1) est bien posé si pour toute donnée f il admet une solution unique u, et cette solution u dépend continûment de la donnée f. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 10 / 82  Introduction Classification des équations aux dérivées partielles Nous allons brièvement classifier les équations aux dérivées partielles linéaires du deuxième ordre portant sur des fonctions réelles de deux variables réelles u(x , y ). Une telle équation s’écrit : ∂2u ∂2u ∂2u ∂u ∂u a 2 + b + c 2 +d +e + fu = g. (2) ∂x ∂x ∂y ∂y ∂x ∂y Pour simplifier nous supposons que les coefficients a, b, c, d, e, f sont constants. Definition On dit que l’équation (2) est elliptique si b 2 − 4ac < 0, parabolique si b 2 − 4ac = 0, et hyperbolique si b 2 − 4ac > 0. Si on applique cette définition aux divers modèles du deuxième ordre que nous avons évoqué dans cette introduction (en remplaçant le couple (x , y ) par les variables (t, x ) en une dimension d’espace), nous concluons que l’équation de la chaleur est parabolique, que l’équation des ondes est hyperbolique, et que l’équation de Laplace est elliptique. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 11 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Principe de la méthode en dimension un On considère le problème unidimensionnel −u 00 (x ) = f (x ), ∀x ∈]0, 1[ ( (3) u(0) = u(1) = 0, où f ∈ C ([0, 1]). Les conditions aux limites considérées ici sont dites de type Dirichlet homogène (le terme homogène désigne les conditions nulles). Cette équation modélise par exemple la diffusion de la chaleur dans un barreau conducteur chauffé (terme source f ) dont les deux extrémités sont plongées dans de la glace. Méthode de différences finies : On cherche une valeur approchée de u en certains points xi de l’intervalle [0, 1]. Ces points sont appelés points de discrétisation du maillage ; on les choisit (pour simplifier) équirépartis, i.e. de la forme : xi = ih, i ∈ {0, · · · , N + 1}, où h = 1/(N + 1) est le pas du maillage (destiné à tendre vers 0 quand le nombre N + 2 de points de discrétisation devient très grand). Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 12 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies On a en particulier aux extrémités : x0 = 0 et xn+1 = 1 ; les autres points xi pour i ∈ {1, · · · , N}, sont dits internes. Figure – Les points de discrétisation du maillage Le principe de la méthode des différences finies consiste à écrire l’équation aux dérivées partielles (3) aux points de discrétisation xi : −u 00 (xi ) = f (xi ) ∀i = 0, · · · , N puis à approcher l’opérateur différentiel (ici −u 00 ) par un quotient différentiel, de manière à en déduire un système d’équations en fonction d’inconnues discrètes censées représenter des approximations de u aux points de discrétisation. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 13 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Voici comment on procède pour l’équation de Poisson unidimensionnelle. Effectuons d’abord un développement de Taylor en xi , en supposant que u ∈ C 4 ([0, 1]) : h2 00 h3 (3) h4 (4) u(xi+1 ) = u(xi ) + hu 0 (xi ) + u (xi ) + u (xi ) + u (ξi ) (4) 2 6 24 h2 h3 (3) h4 (4) u(xi−1 ) = u(xi ) − hu 0 (xi ) + u 00 (xi ) − u (xi ) + u (ηi ) (5) 2 6 24 avec ξi ∈ [xi , xi+1 ] et ηi ∈ [xi−1 , xi ]. En additionnant, on obtient : u(xi+1 ) + u(xi−1 ) = 2u(xi ) + h2 u 00 (xi ) + O(h4 ). Il semble donc raisonnable d’approcher la dérivée seconde −u 00 (xi ) par le quotient différentiel : u(xi+1 ) − 2u(xi ) + u(xi−1 ) −. h2 Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 14 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Sous des hypothèses de régularité sur u, on peut montrer que cette approximation est d’ordre 2 au sens : u(xi+1 ) − 2u(xi ) + u(xi−1 ) Ri = u 00 (xi ) − = O(h2 ). h2 On appelle erreur de consistance au point xi la quantité Ri. L’approximation de u 00 (xi ) par un quotient différentiel suggère de considérer les équations discrètes suivantes : ui+1 − 2ui + ui−1 − = f (xi ) ∀i = 1, · · · , N, h2 dont les inconnues discrètes sont les ui , i = 1, · · · , N. Notons que la première équation fait intervenir u0 tandis que la dernière fait intervenir uN+1. Ces valeurs ne sont pas à proprement parler des inconnues, puisqu’elles sont données par les conditions aux limites. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 15 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies On pose donc u0 = 0 et uN+1 = 0. Le système complet d’équations s’écrit donc ui+1 − 2ui + ui−1  − = f (xi ) ∀i = 1, · · · , N   h2  u0 = 0 (6)    u N+1 = 0. Remarque Attention à ne pas confondre ui et u(xi ) : les équations discrètes (6) font intervenir les inconnues discrètes ui , i = 1, · · · , N et non pas les valeurs u(xi ), i = 1, · · · , N de la solution exacte. En général, ces valeurs ne sont pas les mêmes. Si la discrétisation a été effectuée correctement (comme c’est le cas ici et comme nous le démontrerons mathématiquement plus loin), la résolution du système discret nous permettra d’obtenir des valeurs ui , i = 1, · · · , N des inconnues discrètes qui seront des bonnes approximations des valeurs u(xi ), i = 1, · · · , N de la solution exacte que nous ne pouvons pas, dans le cas général, calculer explicitement. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 16 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Autres conditions aux limites Conditions de Dirichlet non homogènes : Supposons que les conditions aux limites en 0 et en 1 soit maintenant de type Dirichlet non homogènes, c’est-à-dire : u(0) = a et u(1) = b (7) avec a et b pas forcément nuls. Dans ce cas, les équations discrètes du schéma aux différences finies (6) restent identiques mais les valeurs u0 et uN+1 sont maintenant données par u0 = a et uN+1 = b. Conditions de Neumann : On appelle condition de Neumann une condition qui impose une valeur de la dérivée, par exemple : u 0 (0) = a. (8) Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 17 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Conditions de Fourier : On appelle condition de Fourier ou condition de Robin une condition qui impose une relation entre la valeur de la dérivée et la valeur de la solution, par exemple, u 0 (1) + αu(1) = b, (9) avec α > 0. Cette condition est donc un mélange des conditions de Dirichlet et de Neumann et est souvent utilisée pour exprimer une condition de transfert, thermique par exemple, entre un milieu et l’extérieur. Conditions mixtes : on dit que les conditions aux limites sont mixtes si elles sont de type différent sur des portions de frontière du domaine : on a des conditions mixtes dans le cas unidimensionnel si, par exemple, on a une condition de Dirichlet en 0 et une condition de Neumann en 1. Exemple : Prenons par exemple le cas de conditions mixtes, en considérant l’équation (3) en x = 0 avec les conditions (8) et (9) en x = 1. Voyons comment tenir compte de ces nouvelles conditions limites : Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 18 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Pour approcher la condition de Neumann en 0, on effectue un développement de Taylor à l’ordre 1 en 0 : u(x1 ) − u(0) u 0 (0) = + ε(h) = a. h Ceci suggère l’équation discrète suivante pour u0 en écrivant que : u1 − u0 = a ↔ u0 = u1 − ah. h De la même manière, on écrit un développement limité pour la dérivée dans la condition de Fourier (9) : u(1) − u(xN ) + ε(h) + αu(1) = b h ce qui suggère l’approximation suivante : uN+1 − uN uN + bh + αuN+1 = b ↔ uN+1 =. h 1 + αh Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 19 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Principe de la méthode en dimension deux On cherche à approcher la solution u du problème suivant : ( −∆u = f , sur Ω =]0, 1[×]0, 1[, u = 0, sur ∂Ω, où f est une fonction donnée, continue sur [0, 1]. Pour cela, on recouvre Ω de petits rectangles élémentaires de taille h1 = 1/(N1 + 1) dans la direction x1 et h2 = 1/(N2 + 1) dans la direction x2. Figure – Maillage 2d Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 20 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies On cherche, pour chaque indice i1 ∈ {0, · · · , N1 + 1}, i2 ∈ {0, · · · , N2 + 1}, une approximation, notée ui1 ,i2 , de u(i1 h1 , i2 h2 ). En effectuant les développements limités de Taylor dans les deux directions, ∂2u ∂2u on approche 2 (x1 , ·) (respectivement (·, x2 )) par : ∂x1 ∂x22 u(x1 + h1 , ·) − 2u(x1 , ·) + u(x1 − h1 , ·) h12 u(·, x2 + h2 ) − 2u(·, x2 ) + u(·, x2 − h2 ) (respectivement par ). h22 ∂2u ∂2u Compte tenu du fait que ∆u = + , cela donne le schéma suivant : ∂x12 ∂x22 i1 +1,i2 − 2ui1 ,i2 + ui1 −1,i2 ui ,i +1 − 2ui1 ,i2 + ui1 ,i2 −1  u  −   2 − 1 2 2 h1 h2   = f (xi1 ,i2 ) , pour i1 ∈ {1,... , N1 } , i2 ∈ {1,... , N2 } , ui1 ,i2 = 0, pour i1 ∈ {0, N1 + 1} ou i2 ∈ {0, N2 + 1}.  Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 21 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Ce schéma s’appelle usuellement «schéma à cinq points du laplacien» ; il s’agit d’un schéma centré, où pour évaluer la valeur de u au point xi1 ,i2 , on utilise les valeurs en cinq points centrés autour du point xi1 ,i2 : le point xi1 ,i2 lui-même et les quatre points suivants : xi1 ,i2 −1 , xi1 ,i2 +1 , xi1 −1,i2 et xi1 +1,i2. Figure – Le schéma à cinq points du laplacien Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 22 / 82  Equations elliptiques en dimension un et deux Principe de la méthode de différences finies Questions d’analyse numérique Voici un certain nombre de questions du domaine de l’analyse numérique, auxquelles nous tenterons de répondre dans la suite : 1 Le problème obtenu avec des inconnues localisées aux nœuds du maillage admet-il une (unique) solution ? 2 La solution du problème discret satisfait-elle les propriétés physiques qui sont vérifiées par la solution du modèle mathématique ? 3 La solution du problème discret converge-t-elle vers la solution du problème continu lorsque le pas du maillage h tend vers 0 ? Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 23 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies On cherche à discrétiser le problème aux limites suivant : −u 00 (x ) + c(x )u(x ) = f (x ), ( ∀x ∈]0, 1[ (10) u(0) = u(1) = 0, où c ∈ C ([0, 1], R+ ) et f ∈ C ([0, 1], R), qui peut modéliser par exemple un phénomène de diffusion-réaction d’une espèce chimique. On se donne un pas du maillage constant h = 1/(N + 1) et une subdivision de [0, 1], notée (xk )k=0,··· ,N+1 avec x0 = 0 < x1 < · · · < xN < xN+1 = 1. Soit ui l’inconnue discrète associée au nœud i, i = 1, · · · , N. On pose u0 = uN+1 = 0. On obtient les équations discrètes en approchant u 00 (xi ) par quotient différentiel. Cela donne le système suivant : ui+1 − 2ui + ui−1  − + ci ui = fi ∀i = 1, · · · , N   h2  u0 = 0 (11)    u N+1 = 0, avec ci = c(xi ) et fi = f (xi ). Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 24 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Remarque On dit usuellement qu’on a discrétisé le problème par une méthode de différences finies utilisant le «schéma à trois points» de la dérivée seconde. Il s’agit d’un schéma dit «centré» : pour évaluer la valeur de u au point xi , on utilise les valeurs en trois points centrés autour de xi , le point xi lui-même et les points xi+1 et xi−1 , comme schématisé sur cette figure : Figure – Le schéma à 3 points Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 25 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies On peut écrire ces équations sous forme matricielle : Ah Uh = Fh (12) avec 2 + c1 h2   −1 0 ··· 0 ....    −1 2 + c2 h2 −1..   1 ......  Ah = 2  0.. 0. , h   .... .. −1 2 + cN−1 h2 −1     0 ··· 0 −1 2 + cN h2 u1 f1     .. ..  Uh = .  et Fh =  .  uN fN Les questions suivantes surgissent alors naturellement : 1 Le système (12) admet-il une unique solution ? 2 A-t-on convergence de U vers u et en quel sens ? h Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 26 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Proposition Soit c = (c1 , · · · , cN ) ∈ RN tel que ci ≥ 0 pour i = 1, · · · , N, alors la matrice Ah définie dans (12) est symétrique définie positive et donc inversible. Remarque Puisque Ah est symétrique définie positive, donc inversible, cela entraîne l’existence et l’unicité de la solution de (12). Introduisons à présent une autre notion utile : celle de matrice monotone. Definition (Matrice monotone) Soit A ∈ Mn (R) de coefficients ai,j , i, j = 1, · · · , N. On dit que A est à coefficients positifs (ou A ≥ 0) si ai,j ≥ 0, ∀i, j = 1, · · · , N. On dit que A est monotone si A est inversible et A−1 ≥ 0, c-à-d que tous les coefficients de A−1 sont positifs ou nuls. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 27 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies On a la caractérisation suivante des matrices monotones : Proposition Soit A ∈ Mn (R) de coefficients ai,j , i, j = 1, · · · , N. Alors A est une matrice monotone si et seulement si ∀v ∈ RN , Av ≥ 0 ⇒ v ≥ 0. (13) (Les inégalités s’entendent composante par composante.) On introduit une propriété qualitative des solutions du problème (10). Théorème (Principe du maximum) Supposons que c ≥ 0 et que le problème (10) admet une solution u ∈ C 2. Si f ≥ 0 dans ]0, 1[, alors on a : u ≥ 0 dans [0, 1]. Explication : Cette propriété mathématique correspond à l’intuition physique : par exemple, si on chauffe f ≥ 0, la température intérieure sera supérieure à celle au bord u ≥ 0. Il est donc souhaitable que la solution approchée satisfasse la même propriété. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 28 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Lemme (Monotonie de la matrice de discrétisation) Soit = (c1 , · · · , cN ) ∈ RN et Ah ∈ MN (R) définie en (12). Si ci ≥ 0, ∀i = 1, · · · , N alors Ah est une matrice monotone. Explication : L’avantage des schémas dont la matrice est monotone est de satisfaire la propriété de conservation de la positivité qui peut être cruciale dans les applications physiques : Si on a : Ah Uh = Fh avec Fh ≥ 0, alors Uh ≥ 0. D La troisième question concerne la convergence de la méthode : en d’autres termes, a-t-on par exemple, en chacun des sommets du maillage : ui − u(xi ) → 0 quand le pas de discrétisation h tend vers zéro ? Peut on évaluer par ailleurs la précision de la méthode ? Une première étape consiste à évaluer l’erreur dite de «consistance» de ce schéma. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 29 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Definition (Erreur de consistance) On appelle erreur de consistance la quantité obtenue en remplaçant l’inconnue par la solution exacte dans le schéma numérique. Dans le cas du schéma (11), l’erreur de consistance au point xi est donc définie par : u(xi+1 ) − 2u(xi ) + u(xi−1 ) Ri = − + c(xi )u(xi ) − f (xi ). (14) h2 L’erreur de consistance Ri est donc l’erreur qu’on commet en remplaçant l’opérateur −u 00 par le quotient différentiel u(xi+1 ) − 2u(xi ) + u(xi−1 ) −. h2 Cette erreur peut être évaluée si u est suffisamment régulière en effectuant des développements de Taylor. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 30 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Definition (Ordre du schéma) On dit qu’un schéma de discrétisation à N points est d’ordre p s’il existe C ∈ R ne dépendant que de la solution exacte, tel que l’erreur de consistance satisfait : max (Ri ) ≤ Chp , i=1,··· ,N où h est le le pas du maillage. Definition (Consistance du schéma) On dit qu’un schéma de discrétisation est consistant si max (Ri ) → 0 lorsque h → 0 i=1,··· ,N où N est le nombre de points de discrétisation. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 31 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Lemme Si la solution de (10) vérifie u ∈ C 4 ([0, 1]) (ce qui suppose f de classe C 2 ), le schéma (11) est consistant d’ordre 2 et on a plus précisément : h2 |Ri | ≤ sup |u (4) | ∀i = 1, · · · , N. (15) 12 [0,1] Remarque Si on note Ūh =: (u(xi ))i=1,··· ,N le vecteur dont les composantes sont les valeurs exactes de la solution de (10) et Uh = (u1 , · · · , uN ) la solution de (12), on a : Rh = Ah Ūh − Fh = Ah (Ūh − Uh ) (16) où Rh ∈ RN est le vecteur de composantes Ri , i = 1, · · · , N, erreur de consistance en xi définie en (14). Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 32 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Remarque Remarquons que, si u est de dérivée quatrième nulle, alors l’erreur de consistance Rh est nulle, de sorte qu’on a : Ah (Uh − Ūh ) = 0. Comme la matrice Ah est inversible, il s’ensuit que Uh = Ūh , ce qui signifie que pour tout i ∈ {0, · · · , N + 1}, nous avons ui = u(xi ). La solution discrète coïncide donc avec la solution exacte en chacun des sommets (internes ou non) du maillage ! Ceci est très spécifique au cas où la solution exacte u se trouve être une fonction polynomiale de degré au plus 3. Bien sûr, en général ui 6= u(xi ). La preuve de convergence du schéma utilise la notion de consistance, ainsi qu’une notion de stabilité, que nous introduisons maintenant : Definition (Stabilité du schéma) On dit que le schéma (11) est stable, au sens où la norme infinie de la solution approchée est bornée par un nombre ne dépendant que de f. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 33 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Proposition La matrice de discrétisation Ah satisfait : 1 kA−1 h k∞ ≤ , (17) 8 inégalité qui peut aussi s’écrire comme une estimation sur les solutions du système (12) : 1 kUh k∞ ≤ kf k∞. (18) 8 Remarque L’inégalité (18) donne une estimation sur les solutions approchées indépendantes du pas de maillage. C’est ce type d’estimation que l’on recherchera comme garant de la stabilité d’un schéma numérique. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 34 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Nous traçons sur la figure ci-dessous la solution exacte ū(x ) = x (1 − x )/2 et la solution approché Uh dans le cas particulier utilisé dans la preuve de l’estimation de stabilité lorsque la solution discrète coïncide avec la solution exacte aux points de discrétisations. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 35 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Definition (Erreur de discrétisation) On appelle erreur de discrétisation en xi , la différence entre la solution exacte en xi et la i-ème composante de la solution donnée par le schéma numérique : ei = u(xi ) − ui ∀i = 1, · · · , N. (19) Théorème Soit u la solution exacte de (10). On suppose u ∈ C 4 ([0, 1]). Soit Uh la solution de (11). Alors l’erreur de discrétisation définie par (19) satisfait : h2 (4) max |ei | ≤ ku k∞. i=1,··· ,N 96 Le schéma est donc convergent d’ordre 2. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 36 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Illustrons le résultat précédent avec un exemple numérique. On prend c(x ) = 1000 sin(10πx )2 et f (x ) = 1, et nous calculons les approximations par DF pour N = 19, 29, 49, 99 et 199, c-à-d h = 0, 05, 0, 033, 0, 02, 0, 01 et 0, 005. Notons que la solution exacte, a été calculé avec la méthode des DF avec N = 8000, puisqu’aucune formule explicite semble être disponible. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 37 / 82  Equations elliptiques en dimension un et deux Analyse de la méthode des différences finies Remarque Supposons que l’on cherche à calculer la solution approchée de −u 00 = f. Le second membre f est donc une donnée du problème. Supposons que des erreurs soient commises sur cette donnée (par exemple des erreurs d’arrondi, ou des erreurs de mesure). On obtient alors un nouveau système, qui s’écrit Ah Ũh = Fh + εh où εh représente la discrétisation des erreurs commises sur le second membre. Si on résout Ah Ũh = Fh + εh au lieu de Ah Uh = Fh , l’erreur commise sur la solution du système s’écrit : Eh = Ũh − Uh = A−1 h εh. On en déduit que : 1 kEh k∞ ≤ kεh k∞. 8 On a donc une borne d’erreur sur l’erreur qu’on obtient sur la solution du système par rapport à l’erreur commise sur le second membre. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 37 / 82  Méthode des différences finies pour l’équation de la chaleur On se propose de résoudre numériquement le problème unidimensionnel : trouver u : [0, 1] × [0, T ] → R solution de : ∂2u  ∂u (x , t) − (x , t) = 0 pour (x , t) ∈]0, 1[×]0, T [,   ∂x 2   ∂t u(0, t) = u(1, t) = 0 pour t ∈]0, T [, (20)    u(x , t = 0) = u0 (x ) pour x ∈]0, 1[,  où u(x , t) représente la température au point x et au temps t. Il s’agit de l’équation de la chaleur en une dimension d’espace. Dans toute le suite, nous supposons que le problème (20) admet une solution u suffisamment régulière, et, pour simplifier, que la condition initiale u0 est compatible avec la condition à la limite, c’est-à-dire : u0 (0) = u0 (1) = 0. Le but est de calculer une solution approchée en un certain nombre de points (xi , tn ) du domaine espace-temps, c’est-à-dire de l’ensemble des points (x , t) ∈ [0, 1] × [0, T ]. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 38 / 82 Méthode des différences finies pour l’équation de la chaleur Pour cela, on se donne un ensemble de points tn , n = 0, · · · , M + 1 de l’intervalle [0, T ], et un ensemble de points xi , i = 0, · · · , N + 1 de l’intervalle [0, 1]. Pour simplifier, on considère un pas constant en temps et en espace. Soit h = 1/(N + 1) le pas de discrétisation en espace, et k = T /(M + 1), le pas de discrétisation en temps. On pose alors tn = nk pour n = 0, · · · , M + 1 et xi = ih pour i = 0, · · · , N + 1. On cherche alors en chacun des points (xi , tn ) un nombre, noté uin qui soit une «bonne» approximation de u(xi , tn ). Comme la condition initiale u0 est compatible avec la condition à la limite, on a u00 = u0 (0) = u0 (1) = uN+1 0 = 0. Il est clair qu’aux extrémités de l’intervalle en espace [0, 1], il suffit de prendre la condition à la limite exacte, i.e. d’imposer : u0n = uN+1 n =0 pour tout n ∈ {0, · · · , M + 1}. De cette sorte, à tout instant tn , les seules inconnues du problème sont les valeurs correspondant aux sommets internes du maillage en espace (i.e. les points xi correspondant à i ∈ {1, · · · , N}). Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 39 / 82 Méthode des différences finies pour l’équation de la chaleur De même, à l’instant initial t0 = 0, il suffit de prendre pour solution approchée la valeur de la donnée initiale en chacun des points xi , i.e. d’imposer : ui0 = u0 (xi ), pour tout i ∈ {1, · · · , N}. Les seules valeurs qui restent inconnues à ce stade sont donc uin pour i = 1, · · · , N et n = 1, · · · , M + 1. Ainsi, on cherche à calculer une solution approchée du problème (20) en calculant des valeurs uin , i = 1, · · · , N et n = 1, · · · , M + 1, censées être des approximations des valeurs u(xi , tn ), i = 1, · · · , N, et n = 1, · · · , M + 1. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 40 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps L’approximation en temps par la méthode d’Euler explicite consiste à écrire la première équation de (20) en chaque point xi et temps tn , à ∂u approcher la dérivée en temps (xi , tn ) par le quotient différentiel : ∂t u(xi , tn+1 ) − u(xi , tn ) k ∂2u et la dérivée en espace − (xi , tn ) par le quotient différentiel : ∂x 2 u(xi+1 , tn ) − 2u(xi , tn ) + u(xi−1 , tn ) −. h2 On obtient alors le schéma suivant : uin+1 − uin ui+1n − 2u n + u n i i−1 − 2 = 0, i = 1, · · · N, n = 0, · · · , M, 0 k h ui = u0 (xi ) , i = 1, · · · , N, u0n = uN+1 n = 0, n = 0, · · · , M + 1. (21) Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 41 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps On notera Un le vecteur de RN formé par les N inconnues uin , pour i ∈ {1, · · · , N}. Ainsi les N premières équations du schéma peuvent s’écrire vectoriellement sous la forme suivante : Un+1 − Un + Ah Un = 0, pour n = 0,... , M, k où on a posé   2 −1 0 ··· 0   −1 2 −1 · · · 0   1 ..........  Ah = 2 ..... . (22) h   0 · · · −1 2 −1     0 ··· 0 −1 2 La donnée initiale discrète est également un vecteur :   u0 (x1 )  u0 (x2 )   ∈ RN.   U0 = ...     u0 (xN ) Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 42 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Avec ces notations, le schéma est équivalent à : ( Un+1 = (I − kAh ) Un pour n = 0,... , M, (23) U0 = U0 , où I est la matrice identité d’ordre N. Connaissant U0 , on peut ainsi calculer de proche en proche chacun des Un. Le schéma est dit explicite car on peut ainsi très facilement calculer Un+1 à partir de Un sans avoir à inverser de système linéaire. k On peut également noter l’apparition du facteur 2 qui jouera un rôle h important dans la suite. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 43 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Consistance Definition (Erreur de consistance, Euler explicite) On pose ūin = u(xi , tn ) la valeur exacte de la solution en xi et tn du problème (20). L’erreur de consistance du schéma (21) en (xi , tn ), notée Rin , est définie comme la somme des erreurs de consistance en temps et en espace : Rin = Ren + R i b n , avec i en = ūin+1 − ūin ∂u Ri − (xi , tn ) et (24) k ∂t n n n 2 b n = − ūi+1 − 2ūi + ūi−1 + ∂ u (xi , tn ). R i h2 ∂x 2 Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 44 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Proposition (Consistance, Euler explicite) Le schéma (21) est consistant d’ordre 1 en temps et d’ordre 2 en espace, c’est-à-dire qu’il existe C > 0 ne dépendant que de u tel que l’erreur de consistance Rn = (R1n , · · · , RN n t ) définie par (24) vérifie : kRn k∞ = max |Rin | ≤ C (k + h2 ). i=1,··· ,N Démonstration. On a vu lors de l’étude des problèmes elliptiques que l’erreur de consistance en espace Rb n est d’ordre 2. Un développement de Taylor en i temps donne facilement que R e n est d’ordre 1 en temps. i Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 45 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Quelques remarques On aurait pu approcher différemment la dérivée en temps, par exemple en écrivant : ∂u u(xi , tn ) − u(xi , tn−1 ) (xi , tn ) ≈. ∂t k On obtient alors le schéma suivant : uin − uin−1 ui+1 n − 2u n + u n i i−1 − = 0, i = 1, · · · N, n = 1, · · · , M + 1. k h2 Vectoriellement, ce schéma peut s’écrire sous la forme suivante : Un − Un−1 + Ah Un = 0, pour n = 1,... , M + 1. k ou de manière équivalente : Un+1 − Un + Ah Un+1 = 0, pour n = 0,... , M. (25) k Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 46 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Quelques remarques Le schéma (25) est dit implicite, car la formule ci-dessus n’est pas une simple relation de récurrence. En effet, U n+1 apparaît comme la solution d’une équation une fois que U n est connue. Comme la matrice (I + kAh ) est symétrique définie positive (et donc inversible), alors cette équation admet une unique solution et on obtient : ( Un+1 = (I + kAh )−1 Un pour n = 0,... , M, U0 = U0. La mise en œuvre de la méthode d’Euler implicite nécessite la résolution d’un système linéaire à chaque pas de temps, alors que la méthode explicite est simplement un produit matrice-vecteur à chaque pas de temps. La méthode d’Euler implicite est donc plus coûteuse en calcul que la méthode explicite, mais elle a d’autres avantages comme nous le verrons plus tard. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 47 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Quelques remarques Par des développements de Taylor, on montre de manière analogue que le schéma implicite est consistant d’ordre deux en espace et un en temps. Question : Peut on améliorer la précision en temps ? Réponse : Oui, en centrant par exemple le schéma de discrétisation en temps, c’est-à-dire en écrivant : ∂u u(xi , tn+1 ) − u(xi , tn−1 ) (xi , tn ) ≈. ∂t 2k Le schéma qui en résulte (appelé schéma saute-mouton ou schéma de Richardson) s’écrit vectoriellement : Un+1 − Un−1 + Ah Un = 0, pour n = 1,... , M. 2k Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 48 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Quelques remarques Ce schéma est explicite, puisque Un+1 se calcule directement à partir de Un et Un−1 : ( Un+1 = U n−1 − 2kAh Un pour n = 1,... , M, U0 = U0 , U1 = U1. On note qu’il est nécessaire pour ce schéma de connaître U0 et U1. Le calcul de U1 peut par exemple être fait en utilisant par exemple le schéma explicite (21) pour n = 0 (cette étape est alors seulement précise d’ordre un en temps) ; mais on peut aussi choisir d’autres schémas plus précis, par exemple un schéma de type Runge-Kutta d’ordre deux, ceci afin de ne perdre aucune précision dans cette étape d’initialisation. Remarque Ce schéma est à priori d’ordre 2 en temps. Cependant, il est toujours instable, et donc numériquement inutilisable. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 49 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Stabilité Definition (Stabilité d’un schéma de discrétisation temps-espace) On dit qu’un schéma est L∞ -stable si la solution approchée qu’il donne est bornée en norme infinie, indépendamment du pas du maillage. Proposition (Stabilité, Euler explicite) Si la condition CFL (ou condition de Courant, Friedrichs et Lewy) suivante k 1 λ := 2 ≤ , (26) h 2 est vérifiée, alors le schéma (21) est L∞ -stable au sens où : max |uin | ≤ ku0 k∞. i=1,··· ,N n=1,··· ,M Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 50 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Voici une suite de figures correspondant à la méthode d’Euler explicite k 1 appliqué à u0 (x ) = 4x (1 − x ) dans un cas où 2 = 0, 55 >. Pour h 2 chaque valeur de n, on trace les valeurs uin. L’apparition de l’instabilité se produit entre n = 10 et n = 20 et ne fait qu’empirer. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 51 / 82 Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 52 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Remarque Notons que la condition de stabilité (26) est très contraignante d’un point de vue pratique, car elle impose au pas de temps d’être très petit. Par exemple, pour un pas d’espace de l’ordre de 10−3 , cela impose un pas de temps de l’ordre de 10−6. Ceci engendre trop de calculs : le schéma d’Euler explicite, est donc rarement utilisé pour les problèmes paraboliques dans les codes industriels. D’où l’idée d’essayer de trouver des schémas pour lesquels il soit possible de s’affranchir d’une telle condition. On les appelera schémas «inconditionnellement stables», c’est-à-dire stables sans aucune condition liant les pas de temps et d’espace. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 53 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Convergence Definition (Erreur de discrétisation) Soit u la solution exacte du problème (20) et (uin )n=1,··· ,M+1 i=1,··· ,N la solution du schéma Euler explicite (21). Pour i = 1, · · · , N et n = 1, · · · , M + 1, on note ūin = u(xi , tn ) et on appelle erreur de discrétisation au point (xi , tn ) la quantité ein = ūin − uin. L’erreur de discrétisation associée au schéma (21) est alors ke n k∞ = max |ein |. i=1,··· ,N Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 54 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Théorème (Convergence du schéma d’Euler explicite) Sous la condition de stabilité (26), il existe C > 0 ne dépendant que de u tel que ke n+1 k∞ ≤ ke 0 k∞ + TC (k + h2 ) ∀n = 0, · · · , M. Ainsi, si ke 0 k∞ = 0 alors max |ein | tend vers 0 lorsque k et h tendent i=1,··· ,N vers 0 pour tout n = 1, · · · , M. Le schéma (21) est donc convergent. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 55 / 82  Méthode des différences finies pour l’équation de la chaleur Discrétisation par un schéma d’Euler explicite en temps Nous traçons les valeurs discrètes U M+1 correspondant à T = 1, 1, pour la valeur initiale u0 (x ) = 4x (1 − x ). Chaque courbe correspond à différentes valeurs de N allant de 3 à 30, et choisissant M de telle sorte que k ≈ 0, 49. h2 Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 56 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Rappel Commençons par un petit rappel sur les équations différentielles. On considère le problème de Cauchy : ( y 0 (t) = f (y (t)), ∀t ∈ I y (t0 ) = y0. Soit k un pas (constant) de discrétisation, on rappelle que les schémas d’Euler explicite et implicite pour la discrétisatision de ce problème s’ecrivent : y n+1 − y n Euler explicite : = f (y n ), n ≥ 0 (27) k y n+1 − y n Euler implicite : = f (y n+1 ), n ≥ 0, (28) k avec y 0 = y0. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 57 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Le θ-schéma On rappelle également que le θ-schéma, où θ est un paramètre de l’intervalle [0, 1] s’écrit : y n+1 = y n + kθf (y n+1 ) + k(1 − θ)f (y n ). Remarque Remarquons que pour θ = 0 on retrouve le schéma explicite (27) et pour θ = 1 le schéma implicite (28). On peut facilement adapter le θ-schéma à la résolution des équations paraboliques. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 58 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Le θ-schéma pour la discrétisation en temps du problème (20), avec une discrétisation par différences finies en espace s’écrit :   n+1 uin+1 uin θ ui+1 − 2uin+1 + ui−1 n+1 n − 2u n + u n
− (1 − θ) ui+1 i i−1 = + , k h2 h2 i = 1,... , N, n = 0, · · · , M, ui0 = u0 (xi ) , i = 1,... , N, u0n = n uN+1 = 0, n = 0, · · · , M + 1. (29) Remarque Si θ = 0, on retrouve le schéma d’Euler explicite ; si θ = 1, celui d’Euler implicite. Dans le cas où θ = 1/2 ce schéma s’appelle schéma de Crank-Nicolson. Notons que dès que θ > 0, le schéma est implicite, au sens où on n’a pas d’expression explicite de uin+1 en fonction des ujn. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 59 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Proposition (Consistance du θ-schéma) Le θ-schéma (29) pour la discrétisation du problème (20) est d’ordre 2 en espace. Il est d’ordre 2 en temps si θ = 1/2, et d’ordre 1 sinon. Convergence du schéma d’Euler implicite : Prenons θ = 1 dans le θ-schéma : on obtient le schéma d’Euler implicite : (1 + 2λ)uin+1 − λui−1 n+1 n+1 − λui+1 = uin avec λ = k/h2. (30) Proposition (Stabilité L∞ pour Euler implicite) Si (uin )i=1,··· ,N est solution du schéma (30), alors : max uin+1 ≤ max uin ≤ max ui0. i=1,··· ,N i=1,··· ,N i=1,··· ,N De même : min uin+1 ≥ min uin ≥ min ui0. i=1,··· ,N i=1,··· ,N i=1,··· ,N Le schéma (30) est donc inconditionnellement L∞ -stable. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 60 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Théorème (Convergence, Euler implicite) Soit e n = (e1n , · · · , eNn ) l’erreur de discrétisation, définie par ein = u(xi , tn ) − uin pour i = 1, · · · , N, avec u la solution exacte du n=1,··· ,M+1 problème (20) et (uin )i=1,··· ,N la solution du schéma d’Euler implicite (30). Alors, il existe C > 0 ne dépendant que de u tel que ke n+1 k∞ ≤ ke 0 k∞ + TC (k + h2 ) n = 0, · · · , M. Ainsi, si ke 0 k∞ = 0 alors le schéma est donc convergent d’ordre 1 en temps et 2 en espace. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 61 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Nous traçons sur la figure ci-dessous la solution du schéma d’Euler implicite correspondant à u0 (x ) = 4x (1 − x ) à l’instant T = 1, 1 et N = M allant de 120 à 240 par incréments de 12. k Dans les calculs ci-dessus, 2 varie d’environ 130 à 270, sans aucun effet h néfaste, bien sûr. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 62 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Pour comparer entre les schémas d’Euler implicite et de Crank-Nicolson, nous traçons sur les mêmes graphes, les résultats obtenues par chaque schéma avec les mêmes paramètres de discrétisation h et k fixés partout, et la solution exacte, pour différentes valeurs de n. On prend u0 (x ) = sin(πx )/2 + sin(2πx ). Dans ce cas, la solution exacte est donnée 2 2 par u(x , t) = 1/2 sin(πx )e −π t + sin(2πx )e −4π t , ce qui permet de faire les comparaisons. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 63 / 82 Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 64 / 82 Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 65 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Remarque 1 Les deux schémas sont stables et l’ordre supérieur (donc une meilleure précision), du schéma de Crank-Nicolson est clairement visible pour cette donnée initiale particulière. 2 On va mettre en évidence un défaut du schéma Crank-Nicolson lorsque la donnée initiale est discontinue et le pas de temps de l’ordre de h. Pour cela, par exemple, on va travailler sur ]0, L[×]0, T [ avec L = 3 et T = 0.1 et on choisit comme donnée initiale : u0 (x ) = 1[ L , L ] (x ). 3 2 Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 66 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson Avec un petit pas de temps type CFL (ici k = 0.5h2 ) du schéma explicite, on a la solution attendu avec régularisation de la solution. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 67 / 82  Méthode des différences finies pour l’équation de la chaleur Schéma implicite et schéma de Crank-Nicolson En revanche, le même code avec un pas de temps raisonablement gros (ici k = 0.5h : un pas de temps de l’ordre de h) fait apparaître (pour le schéma de Crank-Nicolson) les artefacts ci-après hérités de la discontinuité. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 68 / 82  Méthode des différences finies pour l’équation des ondes On se propose de résoudre numériquement par un schéma de différences finies le problème unidimensionnel : trouver u : [0, 1] × [0, T ] → R solution de :  2 ∂ u ∂2u (x , t) − (x , t) = 0 pour (x , t) ∈]0, 1[×]0, T [,    ∂t 2 ∂x 2    u(0, t) = u(1, t) = 0 pour t ∈]0, T [, u(x , t = 0) = u0 (x ), ∂u (x , t = 0) = u1 (x )   pour x ∈]0, 1[,  ∂t (31) u0 et u1 étant les données initiales du problème. Il s’agit de l’équation des ondes en une dimension d’espace. Hypothèse Nous supposons les conditions initiales u0 et u1 compatibles avec la condition à la limite, c’est-à-dire : u0 (0) = u0 (1) = 0 et u1 (0) = u1 (1) = 0. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 69 / 82  Méthode des différences finies pour l’équation des ondes On se donne : N + 1 points de discrétisation en espace (N ∈ N) distants de 1 h= > 0; N +1 M + 1 points de discrétisation en temps (M ∈ N), le pas de temps T étant alors défini par : k = > 0; M +1 puis on pose : xi = ih, pour i ∈ {0, · · · , N+1}, tn = nk, pour n ∈ {0, · · · , M+1}. On cherche en chacun des points (xi , tn ) une valeur approchée, notée uin , de u(xi , tn ). Question Comment ? Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 70 / 82  Méthode des différences finies pour l’équation des ondes Aux extrémités de l’intervalle en espace [0, 1], on prend (comme d’habitude) la condition à la limite exacte, i.e. on pose : u0n = uN+1 n = 0, pour tout n ∈ {0, · · · , M + 1}. (32) A tout instant tn , les inconnues du problème sont alors les valeurs uin correspondant aux sommets internes du maillage en espace (i.e. les points xi avec i ∈ {1, · · · , N}) ; on note Un ∈ RN le vecteur formé par ces N inconnues. De même, à l’instant initial t0 = 0, on impose la condition initiale exacte (supposée compatible avec la condition à la limite), i.e. :   u0 (x1 )  0  u0 (x2 )   U =.. . (33).     u0 (xN ) Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 71 / 82  Méthode des différences finies pour l’équation des ondes Mais comme il s’agit d’un problème du second ordre en temps, il est aussi nécessaire de se donner U1. Par exemple, en écrivant que : ∂u u(x , k) − u(x , 0) u1 (x ) = (x , 0) ≈ , ∂t k ce qui donne la valeur suivante pour U1 :   u1 (x1 ) 1   0 u1 (x2 )   U =U +k.. , (34).     u1 (xN ) mais d’autres choix sont possibles. Par la suite, nous supposons U0 et U1 donnés. Il reste à préciser le calcul de Un pour n > 1. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 72 / 82  Méthode des différences finies pour l’équation des ondes Le schéma explicite naturel ∂2 ∂2 On approche chacune des dérivées secondes et par le schéma ∂t 2 ∂x 2 (usuel) centré à trois points. Cela donne le schéma suivant : uin+1 − 2uin + uin−1 ui+1 n − 2u n + u n i i−1 − = 0, (35) k2 h2 i = 1, · · · N, n = 1, · · · , M, couplé aux conditions aux limites (32) et initiales (33)-(34) discrètes décrites précédemment. Vectoriellement, ce système peut s’écrire comme un schéma à deux niveaux de la forme : Un+1 − 2Un + Un−1 + Ah Un = 0, n = 1, · · · , M, (36) k2 avec U0 , U1 donnés et la matrice Ah est encore définie par (22). On vérifie aisément que ce schéma est explicite, consistant du second ordre en temps et en espace. Par conséquent, sa convergence est uniquement une question de stabilité ! Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 73 / 82  Méthode des différences finies pour l’équation des ondes Le schéma explicite naturel Nous traçcons le résultat du schéma explicite avec les marques + et la k solution exacte en ligne continue à la figure ci-dessous lorsque = 1. h Nous prenons u0 = sin(πx ) et u1 = 0. La solution exacte est donnée par u(x , t) = sin(πx ) cos(πt). Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 74 / 82  Méthode des différences finies pour l’équation des ondes Le schéma explicite naturel Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 75 / 82  Méthode des différences finies pour l’équation des ondes Le schéma explicite naturel Passons maintenant à visualiser ce qui se passe après quelques itérations k lorsque > 1, même avec la condition initiale très régulière utilisée h ci-dessus. Nous remarquons l’apparition de l’instabilité entre n = 15 et n = 25 : Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 76 / 82  Méthode des différences finies pour l’équation des ondes Le schéma explicite naturel Le schèma explicite est donc conditionnellement stable avec une condition de stabilité donnée par : k/h ≤ 1. Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 77 / 82  Méthode des différences finies pour l’équation des ondes Le schéma implicite naturel Le schéma suivant est la version totalement implicite naturel du précédent : n+1 uin+1 − 2uin + uin−1 ui+1 − 2uin+1 + ui−1 n+1 − = 0, (37) k2 h2 i = 1, · · · N, n = 1, · · · , M, avec les conditions aux limites et les données initiales habituelles. Sous forme vectorielle, il s’écrit : 1 2 1 I + Ah Un+1 = 2 Un − 2 Un−1 , n = 1, · · · , M.
2 (38) k k k Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 78 / 82  Méthode des différences finies pour l’équation des ondes Le schéma implicite naturel ∂u Pour approcher (xi , 0), on a utilisé la formule de Taylor : ∂t ∂u u(xi , k) = u(xi , 0) + k (xi , 0) + O(k 2 ), ∂t ce qui a conduit à l’approximation : ui1 = u0 (xi ) + ku1 (xi ). Ce choix correspond alors à une approximation d’ordre 1 en temps de ∂u. ∂t Pour avoir de l’ordre 2 en temps (plus cohérent avec le schéma à l’intérieur) on pousse plus loin le développement : k 2 ∂2u u(xi , k) = u0 (xi ) + ku1 (xi ) + (xi , 0) + O(k 3 ), 2 ∂t 2 ∂2u ∂2u et on utilise l’équation = , d’où : ∂t 2 ∂x 2 k 2 ∂2u u(xi , k) = u0 (xi ) + ku1 (xi ) + (xi , 0) + O(k 3 ). 2 ∂x 2 Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 79 / 82  Méthode des différences finies pour l’équation des ondes Le schéma implicite naturel Ce dernier choix conduit cette fois à l’approximation : 0 − 2u 0 + u 0 k 2 ui+1 i i−1 ui1 = u0 (xi ) + ku1 (xi ) + , 2 h2 qui s’écrit encore :   u1 (x1 ) k2  u1 (x2 )  U1 = (I − Ah )U0 + k   .. . (39) 2  .   u1 (xN ) Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 80 / 82  Notes bibliographiques et remarques Notes bibliographiques et remarques 1 On pourrait se contenter de dire qu’on ne vous a presque rien dit : le champ de l’analyse numérique consacré à l’approximation des équations aux dérivées partielles est très vaste ! 2 Signalons que la méthode des éléments finis est de nos jours probablement la plus répandue pour résoudre numériquement des équations aux dérivées partielles : G. Allaire, Analyse numérique et optimisation : une introduction à la modélisation mathématique et à la simulation numérique, Editions Ecole Polytechnique, (2005). B. Lucquin, Équations aux dérivées partielles et leurs approximations : niveau M1, Ellipses Éd. Marketing (2004). 3 Parmi les autres techniques répandues, mentionnons la méthode des volumes finis : R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, Handbook of numerical analysis, vol. 7, p. 713-1018 (2000). 4 D’autres méthodes sont également possibles (méthodes spectrales...). Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 81 / 82  Notes bibliographiques et remarques Merci pour votre attention Jawad SALHI, FST Errachidia Module M617: Analyse numérique des EDP A.U. 2022/2023 82 / 82