Résolution numérique de l`équation de Poisson et application à l
Contents
1. help de Scilab Question 8 Ecrire un programme d finissant la matrice A sous forme de matrice creuse cela permet d conomiser espace m moire et temps de calcul On pourra utiliser les fonctions speye spzeros l op rateur colon exemple A i ou A 2 N sparse Question 9 Ecrire un programme qui calcule la solution v de 8 en utilisant la m thode de Cholesky Le programme utilise la matrice creuse d finie dans la question pr c dente et les fonctions chfact et chsolve disponibles sous Scilab Les param tres du programme sont le pas de discr tisation h et le second membre b de 8 Le programme affiche en sortie la norme relative du r sidu LH Question 10 Tester votre programme sur l exemple 12 et comparer les performances des m thodes de relaxation et de Cholesky On pourra comparer la pr cision du r sultat le temps de calcul CPU fonction timer l espace m moire utilis fonctions stacksize who 6 3 3 Calcul de la d formation d une membrane Dans cette partie vous pouvez utiliser la m thode de votre choix SOR ou Cholesky pour r soudre l quation de Poisson qui d termine la d formation d une membrane Il s agit de tester la convergence de la m thode des diff rences finies sur un exemple o la solution de 1 est explicite puis d tudier un exemple plus complexe Question 11 On d finit u x y 4x 1 x y 1 y 1 cos 2q7x cos 2qry cos 2qnx cos 2qry 18 q2. jh h 1 N 1 et0 lt i j lt N T On note u une approximation de u x yj obtenue par un sch ma aux diff rences finies que nous allons d crire On fixe uo UnN41 5 Uio Win 1 0 conditions aux limites 2 Lorsque u est de classe C on a du 1 ne Ci Vi 7 u rix1 yj 2u x yj ulzi 1 yi O h utiliser la formule de Taylor comme dans le 1er chapitre du cours De m me du 1 gyz r Vi z ua Yj 1 Zuli yj ulzi Yj 1 O R L quation 1 donne donc en x y E U Eiy1 Yj Uri 95 ulzi yj ulzi V1 Auf yj f i yj Oh 3 On remplace donc 1 2 par le probl me approch 1 Ja Eo Uii Uii Uigjyi Uig dus f leny 1 lt ij lt N 4 Uoj UN 1 j Uio Uin 50 0 lt i 7 lt N 5 On d finit le vecteur v U11 U21 UN1 U12 U22 UN2 ULN U2 N lt 3 UNN 6 contenant les valeurs approch es de u calculer On note de m me fij f i yj et b h fi 2a Es JN1 fiz f22 TED N 2 pary fin LENS EE fan 7 Le probl me 4 5 s crit sous la forme Av b 8 avec Mx2 R La matrice s crit par blocs de taille N S I0 0 I S JAI 0 A 0 a 0 2 0 Es S oI 0 0 I S o J d signe la matrice identit d ordre N et S MN R s crit 4 1 0 0 1 4 1 0 S 0 1 0 0 4 0 0 1 4 On peut v rifier que admet comme valeurs propres krh Ara 4 sin S sin Eu 1 lt k l
3. tant un entier gt 1 Calculer f Au Tracer les graphes des fonctions u et f pour q 1 et q 6 voir help plot3d Question 12 On consid re le cas q 1 de la question pr c dente Pour diff rentes valeurs de h r soudre num riquement le syst me 8 avec second membre b d fini par 7 Tracer le graphe de l erreur relative al en fonction de h on rappelle que d contient les valeurs de la solution exacte u aux noeuds et retrouver num riquement le r sultat de convergence 15 On utilisera une chelle logarithmique pour tracer le graphe de l erreur voir help plot2d Question 13 Dans l exemple des questions 11 et 12 comment le choix de h est il li q si on veut obtenir u avec une bonne pr cision Question 14 On mod lise la pression exerc e par un doigt sur la membrane en un point de coordonn es o Yo par 1 moe i 5 si l s Cowl lt r p 0 05 0 sinon Calculer le d placement u de la membrane sous l effet de la pression des 5 doigts d une main On repr sentera le graphe de u fonction plot3d ainsi que ses valeurs en niveaux de gris comme dans la figure 1 fonctions grayplot Sgrayplot Pr ciser les valeurs des param tres choisis dans vos calculs position des doigts pas de discr tisation h Remarque si n cessaire l espace m moire utilisable pour les calculs peut tre augment avec la fonction stacksize A Rappel de quelques commandes utiles en SCILAB Apr s avoir lanc Sc
4. lt N avec des vecteurs propres correspondant kD ee ui sin krhi sin IThj La matrice A est donc inversible et sym trique d finie positive Pour calculer num riquement la solution du syst me 8 nous allons utiliser la m thode de Cholesky et la m thode SOR m thode de relaxation Cette derni re est int ressante ici car la matrice A est creuse ses coefficients non nuls sont localis s sur 5 diagonales 2 3 Compl ment sur la m thode SOR On pose M N avec M 2 Let N 2 D U les matrices D L U tant les parties diagonale triangulaire inf rieure et triangulaire sup rieure de d finies en cours et w le param tre de relaxation Lorsque w 0 2 pour toute condition initiale zo la suite xx x gt 0 d finie par Mary Nzk b 9 converge vers la solution de 8 lorsque k oo La condition w 0 2 est n cessaire et suffisante pour la convergence de la m thode car A est sym trique d finie positive La question est maintenant de conna tre la valeur optimale de w i e celle donnant la convergence la plus rapide de la m thode Cette valeur de w est celle qui minimise le rayon spectral de la matrice Lu MHN Pour la matrice A particuli re intervenant dans 8 les valeurs propres de Lu sont connues explicitement En effet Luz z quivaut Uz Lz A w 1 Dz 1 Ww soit en notant comme pr c demment z 1 lt i j lt N les composantes du vecteu
5. R solution num rique de l quation de Poisson et application l tude de la d formation d une membrane Travaux pratiques du cours de m thodes num riques ENSIMAG 1 re ann e 4 Avril 2008 R sum Ces travaux pratiques concernent la r solution de l quation de Poisson en deux di mensions sur un domaine rectangulaire avec des conditions aux limites de Dirichlet ho mog nes Le probl me discr tis par la m thode des diff rences finies conduit la r so lution d un syst me lin aire de grande taille La matrice du syst me est la fois creuse et sym trique d finie positive Deux m thodes vues en cours sont donc bien adapt es sa r solution l une it rative m thode de relaxation et l autre directe m thode de Cholesky Dans ce TP vous allez programmer la m thode des diff rences finies et la r so lution du syst me lin aire en utilisant ces deux m thodes et comparer leur efficacit sur cet exemple FIG 1 D formation d une membrane rectangulaire sous l effet de la pression de cinq doigts Le mod le utilis est une quation de Poisson pour le d placement de la membrane La d formation est indiqu e en niveaux de gris le blanc correspondant un d placement presque nul et le noir un d placement maximal Le calcul num rique est effectu sur un maillage de 400 x 400 points 1 Informations pratiques BUT du TP Un l ment central du cours de m thodes
6. h 0 Aw c gt wl lt klc gt 0 ind pendant de h 14 I vl O h 15 Montrer que prendre h dans 11 garantit que d xz O h Question 6 Reprendre le calcul de la question 3 avec h 1 65 1077 kmag 1000 pour diff rentes valeurs de w dans l intervalle 1 7 2 Tracer le graphe du nombre d it rations effectu es en sortie du programme en fonction de w Pour quelle valeur de w r alise t on le moins d it rations Comparer cette valeur la valeur th orique 10 Question 7 On propose de tester deux crit res d arr t en remplacement de 11 Le premier est Ie zrli Trl lt e ou k gt kmar 16 I et le second zk 0 lt T lt e o k gt kmar 17 Gex 2x0 244 ce dernier crit re s obtient sous l hypoth se limk oo te 0 1 Comparer l effi cacit des crit res 11 16 et 17 sur des exemples de votre choix 3 2 M thode de Cholesky Dans cette partie du TP nous allons r soudre le syst me 8 en utilisant la m thode de Cholesky Il ne s agira pas de programmer la m thode elle m me c est dire la factorisation A TT et la r solution de syst mes triangulaires mais d utiliser des fonctions pr program m es effectuant les calculs voir plus bas Le but de ces questions est de vous faire conna tre diff rents outils matrices creuses factorisation de Cholesky et de vous habituer utiliser l
7. ilab vous pouvez tester les commandes suivantes help pour ouvrir l aide en scilab help mot cl pour obtenir la description de la fonction mot cl apropos mot cl pour obtenir la liste des pages d aide contenant mot cl Les commandes clear clc et clf permettent d effacer respectivement les donn es mises en m moire l cran de commandes et les figures Elles doivent tre ex cut es r guli rement pour viter les erreurs et lib rer la m moire A 1 Ex cuter sous Scilab Les commandes Scilab peuvent tre tap es directement en ligne par exemple gt x 1 gt A ones 3 4 gt x A ou bien crites dans un fichier de commmandes sce Dans ce cas 1 Ouvrir un fichier intitul par exemple test sce comportant les instructions suivantes Symbole commentaire Programme test sce clc clf clear A ones 3 4 1 A 2 Sous scilab tapez gt exec test sce On peut galement d finir des fichiers de fonctions nomm s sci Pour cela 1 Ouvrir un fichier intitul par exemple carre sci comportant les instructions suivantes Fonction carre sci function d carre x d Xx x endfunction 2 Sous scilab chargez et compilez le fichier carre sci gt getf carre sci si le fichier est dans le r pertoire courant La fonction carre est maintenant d finie sous scilab gt x 0 1 2 3 4 gt carre x Un calcul trop long peut et
8. mpte rendu sera dactylographi Nous conseillons fortement d utiliser le logiciel TEX qui est un outil extr mement utilis pour la r daction d articles scientifiques Ce compte rendu n exc dera pas 10 pages et ne comportera pas de programmes Les programmes Scilab nous seront envoy s par email la lisibilit du code et la pertinence des commentaires seront pris en compte dans la note du TP Remise du rapport et des programmes Le TP est rendre au plus tard le lundi 11 Mai 2009 17h00 Il faudra d poser votre compte rendu imprim dans le casier pr vu cet effet envoyer un fichier pdf de votre compte rendu ainsi que vos fichiers Scilab l adresse Adrien Magni imag fr avec copie tpmnensimag yahoo fr Conseils et contacts Une s ance de questions r ponses concernant le TP est pr vue durant les cours d amphi des 20 et 21 Avril Plus g n ralement n h sitez pas demander des conseils votre enseignant la fin de chaque cours Par ailleurs pour toute question pr cision sur le TP vous pouvez tout moment vous adresser Adrien Magni imag fr Evitez de faire le TP au dernier moment ce travail demande du temps et lancer trop de calculs simultan ment risque de saturer les machines au mauvais moment 2 Description du probl me 2 1 Equation de Poisson pour la d formation d une membrane On consid re une membrane lastique horizontale de forme rectangulaire On met la mem brane sous ten
9. n note ry b Azp ce vecteur est appel e le r sidu de la m thode it rative Ecrire un programme qui calcule x jusqu ce que rx Te lt e ou k gt kmar 11 Les param tres du programme sont le pas de discr tisation h le second membre b de 8 la tol rance sur la norme relative du r sidu le nombre d it rations maximal kmax et le param tre de relaxation w Le programme affiche en sortie le nombre d it rations r alis es et la norme relative du dernier r sidu Question 3 On se propose de tester le programme pr c dent sur un exemple On choisit h 107 et on fize w avec 10 valeur optimale th orique On d finit u x y 16xy 1 x 1 y 12 et uij u x j On consid re le vecteur v d fini par 6 etb Av Calculer l erreur relative 13 la fin de l ex cution de programme pr c dent lorsque 1074 et kmax 1000 Question 4 Montrer que l erreur relative 13 d cro t exponentiellement avec k indication tracer le graphe de l erreur en fonction de k en chelle logarithmique Donner une approxima tion de la vitesse de convergence Question 5 On suppose que la solution u de 1 est de classe C sur 0 1 x 0 1 et on note di u x j On note d le vecteur d fini de mani re analogue 6 en rempla ant u par contenant les valeurs de la solution exacte u aux points x x On admet le r sultat de stabilit 1 h2 et le r sultat de convergence quand
10. num riques est la r solution de grands syst mes lin aires issus de la mod lisation de probl mes physiques Bien qu il soit crucial d tudier th oriquement les m thodes les plus classiques pour conna tre notamment leurs conditions de convergence il est tout aussi important de savoir les utilliser correctement 1 Le but de ce TP est de vous faire toucher du doigt les difficult s pouvant intervenir dans l im pl mentation d une m thode num rique et dans l analyse des r sultats obtenus Le travail sera effectu en bin me uniquement Ce TP sera r alis avec l aide de SCILAB logiciel librement distribu et largement utilis La r daction d un rapport vous demandera d tre pr cis et clair tout en restant concis Contenu du rapport Vous devez r diger un compte rendu de TP dans lequel vous r pondrez toutes les questions de l nonc expliciterez les m thodes employ es pr senterez et commenterez les r sultats obtenus Il n y a qu un seul rapport rendre par bin me La qualit de la r daction de la synth se de l analyse des r sultats obtenus sont des crit res importants pour la note Notez que ce sujet ne constitue que la base de ce qui vous est demand soyez critique par rapport vos r sultats proposez d autres id es solutions ou tests La derni re page de votre compte rendu devra tre une sorte de manuel d utilisation o vous explique rez comment utiliser vos programmes Le co
11. que trac de courbe subdivise la fen tre du graphe en une matrice m x n de sous fen tres et s lectionne la p i me pour dessiner le graphe courant l l ment i j de la matrice correspond au graphe num ro i 1 m J Pour dessiner dans une fen tre graphique Tapez scf 1 par exemple pour que la fen tre graphique courante devienne la fen tre 1 Si celle ci n existe pas elle est cr e Pour superposer deux courbes plot x tG gtx 1 m me discr tisation Pour exporter une figure Dans la fen tre graphique exporter cliquer sur le menu File puis Export Dans la fen tre qui s ouvre alors 1 choisir l extension du fichier image pour mettre vos figures en latex choisir Post cript ps 2 choisir entre color et Blacks White 3 choisir l orientation Attention par d faut l image est en landscape format pay sage 4 Entrer le nom du fichier image myfig par exemple sans l extension 10
12. r z R 4 oN Zii j H 2Zi 1 j F Ne T i j 1 F LARGES 1 lt i j lt N 20 j ZN 1 7 Zi 0 5 Zi N 1 0 0 lt i j lt N 1 Lorsque est une valeur propre non nulle le calcul de et z se simplifie gr ce au changement de variable z j ACHI 2p En effet on obtient 4 i Wii Wi j Wij Wij A w lwi 1 lt ij lt N Vw 4 Wo j WN 1 j Wio win4 0 OSJ SNI autrement dit 4 1 7z w 1 est valeur propre de A Les valeurs propres de A tant connues explicitement on obtient explicitement VX en r solvant une quation du second degr puis on en d duit Cela permet de d terminer la valeur propre de plus grand module et donc le rayon spectral de Lu On trouve alors que p est minimal lorsque 2 w EER 10 3 Questions 3 1 Impl mentation de la m thode de relaxation SOR Question 1 Ecrire une fonction qui prend comme arguments w le second membre b de 9 un vecteur xx et calcule la solution zg de 9 Ce programme ne doit pas stocker les matrices A M ou N qui deviennent rapidement trop volumineuses lorsque h est petit On utilise la place la formule donnant explicitement les composantes de zk en fonction des arguments w b k et des coefficients non nuls de A il y a au maximum 5 coefficients non nuls sur chaque ligne de Par la suite on fixe xo 0 comme condition initiale de 9 On d signe par la norme sur R Question 2 O
13. re arr t en cliquant sur la fenetre control puis abort et stop A 2 Vecteurs et matrices La fa on la plus simple de d finir une matrice n x m en Scilab est d entrer au clavier la liste de ses l ments A aiin im de sl Op rations l mentaires A tester sur des exemples A B somme A B produit A B produit terme terme A 2 quivalent A A A 2 quivalent A A det A d terminant de A A transpos e de A inv A inverse de A VVUVVVVVYV A 3 Fonctions chantillonn es Une fonction peut tre d finie par rapport une discr tisation de la variable x ainsi x 0 0 1 1 correspond une discr tisation par pas de 0 1 de x 0 x 1 soit 11 valeurs On d finit des fonctions sur cette grille discr te par exemple y sin 2 x pi x x cos pi xx Somme de deux sinuso des 2x2 parabole A 4 Trac de courbes Pour tracer une courbe y f x sur l intervalle a b H E ji nombre de points de discr tisation dx b a n 1 pas de la discr tisation x a dx b x est chantillonn entre a et b avec un pas de dx plot2d x f x A noter on pourra galement utiliser la commande x linspace a b n pour d finir la grille discr te Pour rajouter un titre xtitle Graphe de la fonction sin Pour mettre une l gende x 0 0 1 10 plot24 x exp x leg exp Pour tracer plusieurs graphes dans une fen tre La commande subplot m n p plac e avant cha
14. sion en tirant sur ses bords et on exerce sur elle une force transversale La membrane l quilibre correspond une surface z u x y o u est solution d une quation de Poisson bidimensionnelle Au f x y 0 1 x 0 1 1 2 2 2 se avec Au Du Ty Le second membre f x y est la pression exerc e sur la surface divis e par la tension de la membrane On compl te 1 par les conditions aux limites u 0 y u 1 y u x 0 u x 1 0 2 conditions aux limites de Dirichlet homog nes qui signifient que le bord de la membrane est maintenu fixe Lorsque f 0 on a simplement u 0 et la membrane est horizontale elle correspond alors au carr 0 1 x 0 1 On souhaite ici tudier la d formation de la 2 membrane sous l action des forces transverses On peut montrer que le probl me aux limites 1 2 admet une solution unique lorsque f est continue sur 0 1 x 0 1 et nous allons calculer num riquement cette solution Remarque le mod le 1 est en fait valable pour de petites d formations de la membrane Notez aussi que l quation de Poisson appara t dans beaucoup d autres contextes conduction de la chaleur diffusion d esp ces chimiques lectromagn tisme m canique des fluides 2 2 Approximation par diff rences finies On consid re une discr tisation du domaine 0 1 x 0 1 suivant une grille r guli re form e par des points de coordonn es x y avec x i h y
Download Pdf Manuals
Related Search