Trang 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
Đang tạo bản xem trước...
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
Hiển thị 20 trang đầu. Tải tài liệu để xem đầy đủ.
Chưa có bình luận nào. Hãy là người đầu tiên!
Bình luận (0)
Chưa có bình luận nào. Hãy là người đầu tiên!