Page 1
Introduction a` l’Analyse Num´erique Ernst Hairer Travaux Pratiques en collaboration avec Assyr Abdulle Universit´e de Gen`eve Octobre 2001 Section de math´ematiques Case postale 240 CH-1211 Gen`eve 24
Generating preview...
Introduction a` l’Analyse Num´erique Ernst Hairer Travaux Pratiques en collaboration avec Assyr Abdulle Universit´e de Gen`eve Octobre 2001 Section de math´ematiques Case postale 240 CH-1211 Gen`eve 24
Table des mati`eres I Int´egration Num´erique 1 I.1 Formules de quadrature et leur ordre . . . . . . . . . . . . . . . . . . . . . . . . 1 I.2 Etude de l’erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 I.3 Superconvergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 I.4 Polynˆomes orthogonaux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 I.5 Formules de quadrature de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . 13 I.6 Un programme adaptatif – TEGRAL . . . . . . . . . . . . . . . . . . . . . . . 16 I.7 L’epsilon-algorithme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 I.8 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 II Interpolation et Approximation 25 II.1 Diff´erences divis´ees et formule de Newton . . . . . . . . . . . . . . . . . . . . 25 II.2 Erreur de l’interpolation et polynˆomes de Chebyshev . . . . . . . . . . . . . . . 28 II.3 Convergence de l’interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 II.4 Influence des erreurs d’arrondi sur l’interpolation . . . . . . . . . . . . . . . . . 36 II.5 Transformation de Fourier discr`ete et interpolation trigonom´etrique . . . . . . . 39 II.6 Transformation de Fourier rapide (FFT) . . . . . . . . . . . . . . . . . . . . . . 44 III II.7II.8II.9Equations Diff´Interpolation par fonctions splineL’erreur du spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .erentielles Ordinaires 4. . . . . . . . . . . . . . . . . . . . . . . . . 46515457 III.1 Un exemple: probl`eme restreint `a trois corps . . . . . . . . . . . . . . . . . . . 58 III.2 M´ethodes de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 III.3 Construction de m´ethodes d’ordre . . . . . . . . . . . . . . . . . . . . . . . . 61 III.4 Un programme `a pas variables . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 III.5 Convergence des m´ethodes de Runge-Kutta . . . . . . . . . . . . . . . . . . . . 66 III.6 M´ethodes multipas (multistep methods) . . . . . . . . . . . . . . . . . . . . . . 69 III.7 Etude de l’erreur locale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .´ 74 III.8 Stabilit´e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 III.9 Convergence des m´ethodes multipas . . . . . . . . . . . . . . . . . . . . . . . . 77 III.10 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 IV Syst`emes d’Equations Lin´eaires 83 IV.1 Elimination de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 IV.2 Le choix du pivot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 IV.3 La condition d’une matrice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 IV.4 La stabilit´e d’un algorithme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
IV.5 L’algorithme de Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 IV.6 Syst`emes surd´etermin´es – m´ethode des moindres carr´es . . . . . . . . . . . . . 96 IV.7 D´ecomposition QR d’une matrice . . . . . . . . . . . . . . . . . . . . . . . . . 97 IV.8 Etude de l’erreur de la m´ethode des moindres carr´es . . . . . . . . . . . . . . . 100 IV.9 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 V Valeurs et Vecteurs Propres 111 V.1 La condition du calcul des valeurs propres . . . . . . . . . . . . . . . . . . . . . 111 V.2 La m´ethode de la puissance . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 V.3 Transformation sous forme de Hessenberg (ou tridiagonale) . . . . . . . . . . . 118 V.4 M´ethode de bissection pour des matrices tridiagonales . . . . . . . . . . . . . . 120 V.5 L’it´eration orthogonale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 V.6 L’ algorithme QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 V.7 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 VI M´ethodes It´eratives – Equations Non Lin´eaires 133 VI.1 M´ethode des approximations successives . . . . . . . . . . . . . . . . . . . . . 133 VI.2 M´ethodes it´eratives pour syst`emes lin´eaires . . . . . . . . . . . . . . . . . . . . 135 VI.3 M´ethode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 VI.4 M´ethode de Gauss-Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 VI.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 VII Travaux Pratiques 147 VII.1 Introduction au FORTRAN 90 . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 VII.2 Int´egration par la r`egle du Trap`eze et Simpson . . . . . . . . . . . . . . . . . . 148 VII.3 Calcul de racines par la m´ethode de la bissection . . . . . . . . . . . . . . . . . 149 VII.4 Wegral: programme d’int´egration d’ordre 8 . . . . . . . . . . . . . . . . . . . . 149 VII.5 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 VII.6 Interpolation et erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 VII.7 R´esolution d’´equations diff´erentielles et calcul de la trajectoire des plan`etes . . . 154 VII.8 D´ecomposition LR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 VII.92D´ecomposition QR et trajectoire d’un ast´ero¨ıde . . . . . . . . . . . . . . . . . . 157 VII.10 FORTRAN 90/95 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 Avant-propos Cette brochure est une revision des polycopi´es distribu´es pendant le cours “Analyse Num´e- rique” ( heures par semaine) donn´e en 2000/2001. Elle est une version corrig´ee des polycopi´es distribu´es en 1991/92. Durant des ann´ees, ce cours a ´et´e donn´e alternativement par Gerhard Wan- ner et par l’auteur. Je remercie Gerhard pour avoir mis `a ma disposition ses notes de cours (en particulier l’exemple de la section VI.4). En ce lieu, j’aimerais remercier Laurent Jay, Luc Rey-Bellet, St´ephane Cirilli, Pierre Leone, Assyr Abdulle, Michael Hauth, mon fils Martin et des nombreux ´etudiants soit pour leur aide dans la pr´eparation des exercices soit pour la correction des erreurs (typographiques et math´ematiques).
Œuvres g´en´erales sur l’analyse num´erique Il y a un grand assortiment de livres qui introduisent le sujet d’analyse num´erique (voir le rayon 65 `a la biblioth`eque de la section de math´ematiques avec plus de 400 livres). En voici quelques ref´erences. Les num´eros entre chrochets (p. ex. [MA 65/403]) vous permettent de trouver le livre facilement `a la biblioth`eque. K.E. Atkinson (1978): An Introduction to Numerical Analysis. John Wiley & Sons. [MA 65/142] W. Boehm & H. Pautzsch (1993): Numerical Methods. AK Peters. [Ma 65/332] G. Dahlquist & A. Bj¨orck (1969): Numeriska Metoder. CWK Gleerup, Bokfoerlag, Lund, Sweden. Traduc- tion anglaise: Numerical Methods. Prentice-Hall, 1974. [MA 65/84] P. Deuflhard & A. Hohmann (1993): Numerische Mathematik I. Eine algorithmisch orientierte Einf¨uhrung. Walter de Gruyter & Co. Traduction anglaise: Numerical analysis. A first course in scientific computation. Walter de Gruyter & Co., 1995. [MA 65/301, MA 65/309] G. Evans (1995): Practical Numerical Analysis. John Wiley & Sons. [MA 65/374] W. Gautschi (1997): Numerical Analysis. An Introduction. Birkh¨auser. [MA 65/393] G. H¨ammerlin & K.-H. Hoffmann (1991): Numerische Mathematik. Zweite Auflage. Springer. [MA 65/303] M.J. Maron (1982): Numerical analysis: a practical approach. Macmillan. W.H. Press, B.P. Flannery, S.A. Teukolsky & W.T. Vetterling (1986): Numerical Recipes (FORTRAN). The Art of Scientific Computing. Cambridge Univ. Press. [MA 65/239] J. Rappaz & M. Picasso (1998): Introduction `a l’analyse num´erique. Presses polytechniques et universitaires romandes. [MA 65/404] A. Quarteroni, R. Sacco & F. Saleri (2000): Numerical Mathematics. Springer. [MA 65/432] M. Schatzman (1991): Analyse num´erique : cours et exercices pour la licence. Paris : InterEditions. [MA E188] H.R. Schwarz (1986): Numerische Mathematik. B.G. Teubner. [MA 65/249] G.W. Stewart (1996): Afternotes on Numerical Analysis. SIAM. [MA 65/389] J. Stoer & R. Bulirsch (1972): Einf¨uhrung in die Numerische Mathematik. Heidelberger Taschenb¨ucher, Springer. Traduction anglaise: Introduction to Numerical Analysis. Springer, 1980. [MA 65/48] Ch. ¨Uberhuber (1995): Computer-Numerik 1 und 2. Springer.
Chapitre I Int´egration Num´eriquefZ:a[ba;f(bx])!dx:IR; Etant donn´e une fonction continue sur un intervalle born´e f(x) (0.1) on cherche `a calculer l’int´egrale (0.2) Dans ce chapitre, nous allons pr´esenter la th´eorie des formules de quadrature (ordre, ´etude de l’er- reur et de la convergence) et nous d´evelopperons les id´ees n´ecessaires `a l’´ecriture d’un programme qui permette de calculer la valeur de (0.2) `a une pr´ecision donn´ee, ceci pour n’importe quelle fonction . Bibliographie sur ce chapitre H. Brass (1977): Quadraturverfahren. Vandenhoeck & Ruprecht. [MA 65/124] P.J. Davis & P. Rabinowitz (1975): Methods of Numerical Integration. Academic Press, New York. H. Engels (1980): Numerical Quadrature and Cubature. Academic Press. G. Evans (1993): Practical Numerical Integration. John Wiley & Sons. [MA 65/336] V.I. Krylov (1959): Priblizhennoe Vychislenie Integralov. Goz. Izd. Fiz.-Mat. Lit., Moscow. Tra- duction anglaise: Approximate calculation of integrals. Macmillan, 1962. [MA 65/185] R. Piessens, E. de Doncker-Kapenga, C.W. ¨Uberhuber & D.K. Kahaner (1983): QUADPACK. A Subroutine Package for Automatic Integration. Springer Series in Comput. Math.,[a;b] vol. 1. [MA 65/210]a=x<x<x<:::<x=b I.1A.H. Stroud (1974):Formules de quadrature et leur ordreger. [MA 65/89]Numerical quadrature and Solution of Ordinary Differential Equations.0 1Zab2f(x)dx=NNjX=01Zxxjj+1f(x)dx: Sprin- La plupart des algorithmes num´eriques proc`edent comme suit: on subdivise en plusieurs sous-intervalles ( ) et on utilise le fait que (1.1) De cette mani`ere, on est ramen´e au calcul de plusieurs int´egrales pour lesquelles la longueur de l’intervalle d’int´egration est relativement petite. Prenons une de ces int´egrales et notons la longueur
2 f(x) Int´egration Num´erique a=xh0j:=xx1j+1xjx2 xjxj+1xN=b g(t):=fF(IGxj. I.1 –+thZxjxjUne division d’un intervalle en sous-intervalles)j+1f(x)dxZ=1hjZ01f(xj+thj)dt: de l’intervalle par . Un changement de variable nous donne alors0g(t)dt: Notons enfin Z1 . Il reste alors `a calculer une approximation de Exemples. 1. LaZ1gformule du point milieu(0t)gdt(t)dt12g(0)g(1+=2)g:(1): (1.2) 2. Laformule du trap`1(0;0g(0))(1=2eze;g(1=2))(1;g(1)) 0g(t)1 2 2 Ces deux formules (point milieu et trap`de degr´(trois pointsi=(3. On obtient laseZ1)01;gg().(ti=)dt(sformule de Simpson,161)))g(0)i+=40g;(1:,:=:2);s+eze) sont exactes sisi l’on passe une parabole (polynˆg1(1)et si l’on approche l’int´: s0repr´1egrale (1.2) par l’aire sous laesente une droite (polynˆ1 ome de degr´s2 e ) par lesome parabole: s7 Z1g(t)dtsXsbig(ci): 0 1 2 1680, Cotes 1711). PourD´efinition 1.14. En g´ci en´eralisant cette id´Une formule de quadrature `, ee (passer un polynˆces formules sont donn´0 a), on obtient lesetages est donn´´ ome de degr´i=1ees dans le tableau I.1.biformules de Newton-Coteseee par par les points ´equidistants(Newton (1.3) Les sont les nœuds de la formule de quadrature et les en sont les poids.
Int´egration Num´serique 1216 1246 16 bi 3 2 ordre2 2889019187TAB288329075. I.1 –38 Formules de Newton-Cotes28812905038poids28890503218 28890757 28819 nomtrap`eze 3 4 41 216 27 272 27 216 41 Simpson 45 46 840840p1840840840840840p NewtonBoole 6 6Z1 s — D´efinition 1.2s7 On dit que l’ordre de la formule de quadrature (1.3) est8 0g(t)pdt=sXi=1big(ci) deggp2, si la formule est exacte1:Weddle pour tous les polynˆomes de degr´e , c.-`a-d., On voit que les formules du point milieu et du trap`Xi=1sbicqi1=1q poureze sont d’ordreqp=1;2;:::;p:. La formule de Newton-q1(1.4) Cotes `Th´geor`aeme 1.31etages a un ordre´;t;::La formule de quadrature (1.3) a un ordre:;tp1c;:::;c(par d´efinition).R01g(t)dtpour si et seulement sipp=P1ssi=1bigg(c(it))=t(1.5) en montrer la suffisance, on utilise le fait qu’un polynˆn´D´eaire deemonstration.b1;:La n´::;becessit´set que l’int´10BBBBe de (1.5) est une cons´c11segralec12::::::c1sainsi que l’expressionequence de (1.4) si l’on pose1CCCC0BBBBome de degr´bb121CCCC=0BBBB1e1=21CCCC:est une combinaison li-sont lin´. Poureaires enlin´eaire pourEn fixant les nœuds. @cs11c(distincts), la condition (1.5) avecs21:::cpss=1As@bsA@1=sArepr´esente un syst`eme Comme la matrice dans (1.6) est inversible (matrice de Vandermonde), la r´... ... ... ... ... esolution de ce syst`q=1;2;3(1.6)eme nous donne une formule de quadrature d’ordre . Si l’on v´erifie les conditions (1.5) pour la formule de Simpson, on fait une observation int´e- ressante. Par d´efinition, il est ´evident que la condition (1.5) est satisfaite pour , mais on
4 11660034++4646q1212=344++16161134==142456=15: Int´egration Num´erique remarque qu’elle satisfait aussi (1.5) pour4 : 3 bi Elle est donc d’ordredegr´sym´etrique.e 2 mais aussi pour des polynˆ. Par cons´0equent, elle n’est pas seulement exacte pour des polynˆomes de degr´c1c2 ec3. Ceci est une propri´c4cmˆieme=c511cs+1et´e g´iben´ierale d’une formule=bs+1iomes de i2m2 FIG. I.2 – Coefficients et nœuds d’une formule de quadrature symetrique2m1 2m1 Th´D´; voir la fig. I.2) a toujours un ordre pair. C.-`emonstration.geor`1(teme 1.4), elle est automatiquement exacte pour les polynˆUne formule de quadrature sym´gChaque polynˆ(t)=2Cm(2tome de degr´1=2)(t2m11=+2)e g2ma-d., si elle est exacte pour les polynˆ1(tetrique (1)peut ˆomes de degr´etre (tec,i1=2)02:5m. omes de degr´1cs+1pour toutite ecrit sous la formeo`´u est de degr´Z0e1(t1=2)2. Il suffit alors de montrer qu’unem1dtR=01(0t:1=2)bi2(mci1dt1=2)2m1+bs+1i(cs+1i1=2)2m1=0 formule sym´etrique est exacte pour . A cause de la sym´etrie de cette fonction, la valeur exacte vaut Pour une formule de quadrature sym´etrique on a . Donc, l’approximation num´I.2 Etude de l’erreurh=(ba)f=erN(erique dexr)=Zbf(x)dx[a;b]NX1hXsbest ´if(egalement z´xj+cih)ero. Afin d’´equidistants (Ensuite, ´cons par une´¸ Prenons une fonctionetudier l’erreur commise en approchant l’int´etudions l’erreur (en ´exp´erience num´erique) et appliquons une formule de quadrature du paragraphe pr´echelle logarithmique), d´efinie sura: , divisons l’intervalle en plusieurs sous-intervallesj=0egrale par une formule de quadrature, commen-i=1f(x) =N(s1)+ecedent.1 (2.1) en fonction de fe (nombre d’´evaluations de la fonction ; on a fe pour Newton-Cotes). Le nombre fe repr´esente une mesure pour le travail (proportionnel au temps de
Int´egration Num´Z03eriquecos(x)exp(sin(x))dx Z0N2cos=(1x;)2dx:;4;8;16;32;::: 5 calcul sur un ordinateur). La fig. I.3 montre les r´esultats (pour ) obtenus par les formules de Newton-Cotes pour les deux int´egrales : et 103 103 fe fe 102 102 101 244erreur 101 s=6 668erreur log10100hF(IG). I.3 –10−3 Le travail fe en fonction de l’erreur pour les formules de Newton-Cotestrap`Newton (ordreSimpson (ordre10eze (ordre−6 101=p−9 ))) p10−12 100 10−3 Boole (ordreWeddle (ordre10−6 (ordre10log−9))10)(err10)−12 En ´etudiant les r´esultats de la fig. I.3, nous constatons que :Z Explication des r´–– pour un travail ´– la pente de chaque droite estfe d´epend lin´esultats de la fig. I.3.Eequivalent, les formules avec un ordre ´(f;eairement du nombre de chiffres exacts, donn´x0;h)==hxZx000+(o`1hEtudions d’abord l’erreur faite sur un sous-intervalle defuf((xx0est l’ordre de la formule);)+dxth)dthXi=1sXib=1sielev´f(bxif0e ont une meilleure pr´(+x0c+ih)ce parih): ecision.; longueur fE(f;x0;h)=qX0hqqx+1!0Z01tqdtXi=1sbicqiff((qx)(0x+0)th)f(x0+cih)(2.2) s´En supposanteries de Taylor (d´suffisamment diff´evelopp´ees autour de=erentiable, on peut remplacerhpp+1!p1), et on obtient ainsi+111Xi=1sbXiscpif(p)(x0)+pO(hp+2et) p+par les1 C=p!p+1i=1bicpi (2.3) (ici, on a bien suppos´e que la formule de quadrature ait l’ordre mais pas l’ordre ). La constante (2.4)
err=NjX=01E(f;xj;h)ChpNjX=01hf(p)(xj)ChChhppZ+1abf(p)(x)dx=Chpf(p1)(b)f(pO1)(h(pa+2)): 6 Int´egration Num´erique s’appelledans (2.3) soit n´C2constante de l’erreur=hegligeable par rapport au termelog10(err). Supposons quelog10(C1)plogsoit petit de mani`10(, alors on obtienth)Constere `+pa ce que le termelog10():errC1hp etCette formule nous permet de mieux comprendre les r´fe1=p , nous avons log10()esultats de la fig. I.3. Commelog10(err) (2.5) fe Ceci montre la d´deEstimation rigoureuse de l’erreur.reur d’une formule de quadrature. Une telle estimation nous permettra de d´f:[.x0;x0+hependance lin´]!IRE(kf;eaire entrex0;h)Notre but suivant est de trouver une estimation exacte de l’er-=hk+1Z0fe1Netk()f(k)(x0+p, et aussi le fait que la pente soith)demontrer des th´k eor`kemesp de convergence et assurera une certaine pr´N() ecision du r´esultat num´erique. Th´Si eor`kNeme 2.1k()=Consid´(1k!)erons une formule de quadrature d’ordrekestXi=1sfois continˆbi(c(ikument diff´1)!)k+1 erentiable, l’erreur (2.2) v´()k+1:=et un entier(()0kerifie1satisfaisant>00 . (2.6) o`u , le noyau de Peano, est donn´f(x0+th)=EkjX(=0f1;x(th0j;!)hj)fe par(j)(x0)+ho`kuZ0t(t(k)1)!k1f(k)(x0+h)sisid ,. D´pemonstration.k Nous introduisons la s´Z0t(t)k1erie de Taylor avec resteg()d=Z01(t)k+1g1 ()d (2.7) et le fait que la partie polynomiale de (2.7) ne donne pas de contribution `dans la formule (2.2) pourE(f;x0;h)=hk+1Z01Z01(t(. En utilisantk)1)!k+1dtXi=1sbi(c(ik1)!)k+1f(k)(a l’erreur (`x0+h)da cause de: 1:), nous obtenons Une ´evaluation de l’int´egrale int´erieure donne le r´esultat. voir le paragraphe III.7 du livre de E. Hairer & G. Wanner (1995), Analysis by Its History. Undergraduate Texts in Mathematics, Springer.
pInt´egration Num´NNNZkkk10(0)((1)N)p(===)k00deriqueN=k11(kk)1211cckiiXksb012ipcipii===11;;C:::6=:::;;cssik=2 7 Th´a)b)et un nombreeor`Neme 2.2 (Propri´i01=(1);:::;satisfaisantspourpet´!es du noyau de Peano)p+pour1si i=1((pour. Alors, on a:Consid´1);si erons une formule de quadrature d’ordre); bi ci c) pour si ( ); e)d) ( est lin´).eaire par morceaux, de pentec1b1 c2b2c3(constante de l’erreur (2.4));b3et avec des sauts de hauteurc4b4 aux points N1()=0 1 <11==22N1(N)2()=((12=2)2=2 111==22 Les noyaux de Peano pour la formule du point milieu sontFNIG1. I.4 –()Le noyau de Peano Nsd’une formule de quadrature2=(5) si si si si (voir la fig. I.5). Pour la formule de Newton-Cotes ( ), ils sont dessin´es dans la fig. I.6. [a;b] fFIG:[. I.5 –a;erbr]=!Noyaux de Peano pour la formule du point milieuZIaRbfk(x)dxNjX=01hjXi=1sbif(xhjj+=cxihj+1j):xj entierGrˆh=ace au r´max. Pour une division arbitraire (´jesultat du th´hj pjerpreor`jkeme pr´hk(bec´edent, on peut facilement estimer l’erreur pour l’intervalleaequidistante ou non;)Z01jNk()jdxmax2[a;b]jf(k)(x)j), notons l’erreur par(2.8) Th´eor`eme 2.3 Soit fois continˆument diff´erentiable et soit l’ordre de la formule de quadrature ´egal `a ( ). Alors, l’erreur (2.8) admet l’estimation (2.9) o`u .
8 N1() N2() N3()Int´egration Num´erique N4() N5() N6() jE(f;x0;h)jhk+1Z1jNk()jjf(k)(x0+h)jd s=5 D´emonstration.FIGLa formule (2.6) donne. I.6 – Noyaux de Peano pour la formule de Newton-Cotes avechk+1Z001jNk()jdx2[maxx0;x0+h]jf(k)(x)j: jerrjNjX=01jE(f;xj;hj)jNjX=01h|hkj{zk+1}hjZ01jNk()jdx|2[maxxxmaxj2;x[a;bj+1]{z]jjff((kk))((xx))j}j Comme l’erreur (2.8) est la somme des erreurs sur les sous-intervalles de la division, on obtientjerrjPhNj=021(hbj=ab)1amaxjf00(x)j; ce qui montre l’assertion (2.9), carExemples. Pour la formule du point milieu, on ajerrjh2(ba)24121.xxmax22[[a;ba;b]]jf00(x)j; pour la formule du trap`eze jerrjs6h=45(ba)288011xmax2[a;b]jf(4)(6)(x)j; pour la formule de Simpsonjerrjh(ba)1935360xmax2[a;b]jf(x)j: pour la formule de Newton-Cotes ( )
Z1jN()jRd01jN=p(Z)1jNds(=)d5=11321N66+(12)16+3236+716=[0;1]1: Int´Le calcul demule de Newton-Cotes (fig. I.6) et on utilise la propri´egration Num´06 erique0c1;:6::pour ces formules n’est pas difficile. Consid´;et´cse (d) du th´). On constate que6!7eor`eme 2.2. Ceci donne904 90ne change pas de signe sur2 904erons par exemple la for-90 (bi1935360;ci)(voir9 I.3 Superconvergenceps (bi;ci)si=1 bici ps formule de l’exercice 1.Question.Si l’on fixe les nœudsun ordre Y a-t-il un choix particulier desZ0. On obtient les poids1M(st)+g(mt)dt=(distincts), il existe une unique formule de quadrature0M(t)=soit par la r´(tpermettant d’avoir un ordre sup´c1):esolution du syst`::(tgcs()t:)eme lin´eaire (1.6), soit par laerieur?m1 ayant Th´eor`eme 3.1 Soit une formule de quadrature d’ordre et soit Alors, l’ordre estfr(t)s1M(t)fg(t)si et seulement sim1f(t)f(t)=M(t)gs(t+)+mr(t)1 (3.1) pour tout polynˆome de degr´e . (3.2) D´diviseremonstration.par Soit et d’´ZecrireXiun polynˆ=1s01bfi(ft()cdti)ome de degr´==sous la formeZXi=1s01bMiM|(=t{z)(egc0(i)}t)gdt(ci+)+Z0Xi1=1. L’id´sr(bti)rdt(ee, due `ci):a Jacobi (1826), est de o`u deg Zf(t)et deg . Alors, l’int´R01M(t)g(tegrale exacte et l’approximation num´)dt=rs(0=t)3 s4erique satisfont Comme la formule de quadrature est exacte pourexacte pour0=01(tc1si et seulement si)(tc2)(tcc33)=dt1==14=143((c(c1c11+++cc2c2)2=+).3=2c+(l’ordre est3+)c311cc+12c=2(2c1:c2+c1cpar hypoth`3+c2c3)21ese), elle estc1c2c3; Exemple 3.2 Pour qu’une formule de quadrature `p a etages ait un ordre´ s , il faut que ce qui est ´equivalent `a p2s: Corollaire 3.3 Si est l’ordre d’une formule de quadrature `a etages, alors´ (3.3)
10 Z01M(t)M(t)dt=g(Zt0)1(tc1)2::s:(tcs)p2dt>2s0Int´+:egration Num´1 erique D´dans (3.2) est nulle pour tout polynˆmemonstration.1M(t) Supposons, par l’absurde, que l’ordre satisfasses hf;gomei=R01f(de degr´t)gmM(t=)(tdt)es2s. Ceci contredit le fait que. Alors, l’int´egrale Pour construire une formule de quadrature d’ordre (si elle existe), il faut trouver un poly- nˆ!ome:(a;b)!de degr´IRe qui satisfasse (3.2) avec . Ceci peut ˆetre r´ealis´e `a l’aide de polynˆomes orthogonaux, car la condition (3.2) signifie que est orthogonal `a tous les polynˆomes de degr´e pour le produit scalaire . I.4Nous nous mettons dans une situation un peu plus g´!(t)Polynˆ=etomes orthogonaux(0;1)Zab!!(t()tj)tj=k!dt(1t)<>10 en´erale. Consid´tk2=(0a;;1b!;)2(t;erons une)::=:1=p1fonction de poidst2 (1;1) qui satisfasse les deux propri´et´es suivantes: Zbpour (4.1) ouDes exemples typiques sontsur . hf;gsur un intervalle born´i=a!(pourt)f(t)g(t)dt:e, ou sur (4.2), Sur l’ensemble des polynˆfomes (`a coefficients r´geels), nous d´() efinissons unhf;gi=0:produit scalaire par (4.3) pLes propri´de ce produit scalaire, on peut parler d’orthogonalit´0;p1;p2;:et´:es du produit scalaire (bilin´: hppkk;(gt)iest orthogonal `==t0k+earit´ae, sym´etrie et positivit´e. On dit quek1e) sont faciles `k1:a v´erifier. A l’aide(4.4) Th´eor`eme 4.1 (existence et unicit´e des polynˆomes orthogonaux) Il existe une suite de polynˆomes telle quek+1k+1ppk+11((tt))==polynˆpour tout polynˆ(0t;ome de degr´k+1)pp0(kt()tome de degr´)=e1k2+1pke1(t) (normalisation)(orthogonalit´e) Ces polynˆomes sont uniques et ils satisfont la formule de r´k+1=hhtppkk;;ppkkii; k2+1=ecurrencehpkhpk1;;ppkki1i: (4.5) o`u les coefficients , sont donn´es par (4.6)
Int´egration Num´PPTkkk(((tt;)))(t)erique((((a;111;;;b1)1)1)) (11=tp!)11((1t)+t2t) ;>1PTPkkk((1)(1);)(1)==11=k+k 11 notationLH(kk()t()t)(1(0;1;1))TAB. I.2 –teet2tPolynˆomes orthogonauxnomLegendre >1pLnormalisation((kt))(0)=1=k+k k pk+1 pk+1(t)=tpk(t)+ChebyshevJacobi,Laguerre,HermiteXkjpj(t)p;k+1(t)—0p0tp;pk1(;t):::;pk D´propri´sons le polynˆemonstration.et´es demand´omep(Orthogonalisation de Gram-Schmidt). Le polynˆkees. Supposons, par r´. La normalisation choisie implique quek ecurrence, qu’on connaisse d´j=0 omejej`a soit un polynˆsatisfait les deuxp0;et construi-p1;:::ome;pk de degr´e , et on a 0=hpk+1;pki=htpk;pki+khpk;pki: k k=k+1 k+1 (4.7) Il reste `mule (4.7) aveccar chaque polynˆa montrer que l’orthogonalit´ome de degr´hdonnetpk;0p=k1hipek+1peut ˆ;pke d´etre ´1ietermine uniquement les=ecrit comme une combinaison lin´htpk;pk1i+pk1k1hpk1;. Le produit scalaire de la for-pk1i:eaire de . Cette relation nous permet de calculerEnsuite, nous prenons le produit scalaire de (4.7) avechtpk;pkp1ii=i0hp2=k;ftph0p;kk1+1;1:i;:p=:i;ikh; on obtient=pk;2pgkhktp+1k=;piik2+1+: ihpik;piiavec1i=hpk;pd´kefini dans (4.6).i; (4.8) scalaire de (4.7) avecNous reformulonset nous ins´erons (4.9) dans (4.8) pour obteniro`uii=comme suit0 hpk|;tp{zipol. de degr´i=}i0. Finalement, nous prenons le produite pk+1(t) (4.9) et nous en d´r´ecurrence (4.5).Cette construction desLes fonctions de poids et les intervalles consid´eduisons que ne nous laisse pas d’autre choix, d’o`. Avec les valeurs deer´es pour quelques polynˆainsi trouv´M(t)u l’unicit´ees, la formule (4.7) devient lae deomes orthogonaux. importants sont pr´esent´es dans le tableau I.2. Souvent, on utilise une autre normalisation (voir derni`ere colonne) que celle du th´eor`eme 4.1. Les polynˆomes du th´eor`eme 4.1 vont jouer le rˆole de dans (3.2). Nous sommes fort int´eress´es `a trouver de tels polynˆomes ayant uniquement des racines r´eelles.
12 pk(t)pkt(1t;):::;tr kpk(t) Int´egration Num´((a;a;bb)) eriquekpk(1t) Th´eor`eme 4.2 Soit un polynˆome de degr´e orthogonal `a tous les polynˆomes de degr´e . Alors, toutes les racines deD´change de signe (voir la fig. I.7).emonstration. Notons paratpsont r´1k(t)eelles, simples et dans l’intervalle ouvertles racines de qui sont r´t2eelles, dansbg(t) .et o`u (tt1):::(ttr)r0=6=kZarb!|><{z(t0k})p|k(t){zg(t)}dt=hpkr;g<i=k0; g(t)= FIG. I.7 – Illustration de la d´emonstration du th´(a;b) eor`eme 4.2 Le but est de montrer . Supposons, par l’absurde, que . Avec le polynˆome de degr´e , on a ne change pas d’o`u la contradiction.pk pk(t)=Ckde signe sur!1(t)dtdkk!(t)(tka)k(bt)k;Ck Th´eor`eme 4.3 (formule de Rodrigues) Les polynˆomes orthogonaux du th´eor`eme 4.1 satisfont si la partie droite de (4.10) est un polynˆhpk;gg(it=)hpCk;kgZiab=!0(t)!ome de degr´1(t)dtdkk:k:e:. La constante1g(t)dt est d´etermin´ee par la(4.10) normalisation deD´d´efini par (4.10), satisfaitemonstration. Soit. ==un polynˆC::k:|dt=dkk(11ome de degr´1). Plusieurs int´:=k:{z:0CkgZ(atb)eba}::egrations par parties donnent:CkgZ(ka)b(. Il suffit de montrer que le polynˆdttd)kkdt11=:0:;:g0(t)dt ome, ce qui montre l’affirmation.
Int´egration Num´Pk(1)erique=1Pkd(kkt)(1+t)k(1t)k=(P(1)k(1tk;)1)k!21k; !(t)=131 PolynˆAu lieu de supposer que le coefficient le plus haut dedemande queomes de Legendret(voir le tableau I.2). Comme=dt1. On consid`1(=1)Cere l’intervallekkd(kt1)=1kk!k2ksoitet la fonction de poidsk(voir le th´eor`eme 4.1), on. la formule (4.10) pourLegendre sont donn´P0(t)=es par1; PP1kdonne((tt))==t;k!2kP2dt(tk)=(132+t2t)(112;. Par cons´t)P:3(t)equent, les polynˆ=25t332t:omes de (4.11) Les premiers de ces polynˆomes sont P3(t)P4(t)PP51((tt)) (4.12) Ils sont dessin´es dans la fig. I.8. P2(t) (1+PPt)kkk(((1tt))==Pt)kkP(k=(t(1)t)t2)kkk FIG. I.8 – Polynˆomes de Legendre Une cons´equence desZ11M(t)g(t)dt=0M(t)sisiest queest pairest impair.gss1: p=2s(4.13) I.5Dans ce paragraphe, nous construisons des formules de quadrature avec un ordrele th´polynˆeor`omes de degr´Formules de quadrature de Gausseme 4.1, il existe un unique polynˆe 0, c.-`a-d.,M(t)=omeCPs(2tde degr´1);e qui soit orthogonal `a tous les. Selon si deg Ce polynˆome est donn´e par (5.1)
14 Ps(x) sZP01s(2Pst(2t1)1)g(t)dt=12Z11Ps(x)g((x+1)=2)dx(0;=1)Int´0:egration Num´erique o`u est le eme` polynˆome de Legendre, car 2s Toutes les racines dele th´cbeor`11;;::eme 3.1 nous donne le r´::::;;bcss esultat suivant.sont r´Ps(2eelles et situ´t1)sees dans l’intervalle (th´eor`eme 4.2). Alors, Th´seor`3eme 5.1 Il existe une formule de quadrature unique d’ordre . Elle est donn´ee par: sont les racines de ; Exemple 5.2 (formules de Gauss)formules de Gaussssssont donn´===123. Pour de petites valeurs de:::es par (1.6).ZZZ00111ggg(((ttt)))dtdtdtLes formules de quadrature du th´g125gg12121, elles sont faciles `p6p315++12g812g+1pa obtenir: il suffit de calculer les6+3eor`5eme pr´g1+ec´pedent s’appellent15 racines de (4.12) et de r´esoudre le syst`0 eme (1.6) tout en exploitant la sym´18210 182 18etrie de la formule. Pour210 , on obtient ainsi: s s(k10+1)P(x)=(2k+1)x(formule du point milieu)P(x)Ps(kxP)(x); Sir´esolution exacte du syst`1;:est grand (disons::;s Peme (1.6) peut poser des probl`s(x), le calcul exact des racines de)k+1 ci=(1+Psemes. D´(kxi))=2ecrivons alors leur calcul pratique.kn’est pas toujours possible et lax1 Calcul de nœuds. En utilisant la formule de r´ecurrence les nœuds de la formule de Gauss `(exercice 11) on peut facilement calculer la valeur dedu polynˆome bipeut alors ˆ=(1a l’aide dei2etre fait par bissection (voir exercice??? 12), et on obtient)(1Ps0(i))2=s2(P1pour uns.1(i2i))2;donn´e. Le calcul des racines(5.2) Calcul de poids. Au lieu de r´esoudre le syst`eme (1.6), on peut aussi utiliser la formule explicitep2sk 1ks l’exercice 11.qui est donn´cee ici sans d´1;:M::(;tc)s=emonstration00;P:s:(2:;tk1)0+26=. La deuxi`01Ps1(2teme identit´1)+b:1:;::e de (5.3) est une cons´+::;bksPsk(2t1) equence de(5.3) Th´si les nœuds2eor`:eme 5.3 Une formule de quadrature a l’ordresont les racines de (avec ) si et seulement(5.4) (avec des nombres r´eels , ) et si les poids sont donn´es par (1.6). voir M. Abramowitz & I.A. Stegun, Handbook of Mathematical Functions, Dover Publ., page 887
Int´egration Num´erique M(t)=0PsM(2(tt)=1)(t+Xi=1sc1)iP:s::si(2(+tt(sc1)s):k) M(15t) D´emonstration. Ecrivons le polynˆ´ ome sc1;:k::;c1s k+1sous la forme=:::=s=0 D’apr`est orthogonal aux polynˆd´emontre (5.4).Pk(1)es le th´=s1eor`=2eme 3.1, la formule de quadrature a l’orderPk(omes de degr´1)=M((1)t)k=ePs(2ct1s2==s1)03. Ceci est le cas si2Pcss=2(2s1t=41); si et seulement si, ce qui Exemple 5.4 (formules de Lobatto)on obtient une formule de quadrature d’ordreZ01g(t)dt121g(0)+Si125g12p10sont les racines de5. Elle s’appelle+125g12+p10formule de Lobatto5+121g(1): . A cause(5.5) de et de , on a et . Des cas particuliers sont la formuleZ du trap`Etudions encore laeze ( ) et la formule de Simpson (g(tpositivit´)0e de l’approximation num´(0;1) =). Pour)erique obtenue par une formule de qua-X0s1, on obtientg(t)dt0: drature. On sait que l’int´g(t)egrale exacte satisfaitg0(t)sur b(0i>;1)0 i==)1;:::;si=1big(ci)0 (5.6) L’analogue num´erique sur 2s1s (5.7) est vrai pour toute fonction(on ne consid`bi ere pas de formule de quadrature ayant des poids nuls).si et seulement si`i(t)=Y((ctisccj1j)) (5.8) poidsD´Th´emonstration.eor`eme 5.5satisfont (5.8).Pour une formule de quadrature d’ordreConsid´erons le polynˆbi=jX=1some de degr´bj`2i(cj)j6==iZe01`2i(t)dt>0(;etant le nombre d’´´ etages), les et int´egrons son carr´e par la formule de quadrature. On obtient`2i(t) car la formule de quadrature est exacte pour .
16 s8 ss=10 Int´egration Num´erique Remarque.le paragraphe I.1) ne satisfont (5.8) que pourI.6(a;b)s=Un programme adaptatif – TEGRAL15Les formules de Gauss satisfont (5.8) pour toutf:[a;b]!IRI=faet pour=x0<. Les formules de Newton-Cotes (voirx1<.:::R<abfx(Nx)=dxbg 30 Posons-nous le probl`qui, pour une fonctionrelative de TOL. Si l’on fixe la formule de quadrature (par exemple la formule de Gauss d’ordrefeme d’´Rbecrire un programmeIFUNCTION TEGRAL ( FCN, A, B, TOL )Zabf(xdonn´)dxee, calcule la valeur de(a;bZ)abjf(x)jdx: a une pr´` ecision avec ), il faut trouver une divisionaf(x)dx de l’intervalle afin que l’approximation num´erique satisfasse relativeennuis dans le cas o`Pour une fonctionest born´ee parune changeant pas de signe surTOL. On a mis la valeur absolue sous l’int´Iest tr`Rabfes petit ou nul.(x)dx TOLs(x, la condition (6.1) signifie que0;x0+hegrale de droite pour ´)(a;b) eviter desl’erreur(6.1) D´etermination de la division.– choix de la division pour que (6.1) soit satisfait;– estimation de l’erreurPour ´ecrire un tel programme, on est confront´Pour un sous-intervalle((xx00;;xx00++hh.))==e aux deux probl`hhXXii=1=1sbbiifjf((xx00++cciemes suivants:ihh));dej: , on sait calculer les valeurs (a;b)(xresabs0;xres(0a;+b)h)(a;(xb0);x0+h)Zxx00+hf(x)dx: (6.2)(6.3) Supposons, pour le moment, qu’on connaisse aussi une estimation de l’erreurj(a;b)j (a;b); L’algorithme pour trouver une division convenable est le suivant:i) on calcule res((Ia;1)(ba;)err, bresabs)(I1X)Net(Ierr1res) . SiI(1I2=)Rab(Xfa;N(x(a)dx+(I2b))=2)(II2)2=((a+bN)=2=;b2)(6.4) ii) on subdivisecalculeon accepteresres , resabsen deux sous-intervallescomme approximation dej=1, errjerr((II1j))etj+res TOL(I2), resabsjresabs=1 et on arrˆ,(errIj)et:ete le calcul; sinon. On poseRabf(x)dxet onet on regarde si err TOL resabs (6.5) Si (6.5) est v´erifi´e, on accepte res res comme approximation de ; sinon
Showing the first 20 pages. Download the document to see it all.
No comments yet. Be the first!
Comments (0)
No comments yet. Be the first!