Elemente de Simulare Numerica v3.1

ELEMENTE DE SIMULARE NUMERICĂ ÎN MATLAB

Curs 1

1

OBIECTIVE

Însușirea principalelor tehnici de simulare nu-merică și a unor metode avansate de calcul numeric al soluțiilor unor clase de probleme ma-tematice.

Folosirea mediului de programare MATLAB pen-tru rezolvarea pe calculator a unor probleme cu grad ridicat de complexitate.

2

Cuprins. Teme abordate

1. Ecuații neliniare: metoda lui Newton, algoritmi pentru aflarea punctelor fixe.

2. Derivarea și integrarea numerică: aproximarea derivatelor funcțiilor, metoda trapezelor și metoda lui Simpson de aproximare a integralei.

3. Sisteme liniare: metode de factorizare și metode iterative.

4. Vectori proprii și valori proprii: metoda puterilor.

5. Aproximarea funcțiilor și datelor: interpolare, functii spline, metoda celor mai mici pătrate.

6. Ecuații diferențiale ordinare: metode de tip Euler, probleme de stabilitate, metode de ordin superior, metode predictor-corector, sisteme de ecuații.

7. Metode numerice pentru probleme la limită: aproximare prin dife-rențe finite, metoda elementelor finite pentru dimensiunile 1 și 2.

3

Bibliografie

• A. Quarteroni, F. Saleri — Scientific Computing with MATLAB and Octave, 2nd Edition, Springer-Verlag, Berlin, 2006.

• A. Quarteroni, R. Sacco, F. Saleri — Numerical Mathematics, Springer-Verlag, New York, 2000.

• M. Ghinea, V. Fireteanu — MATLAB. Calcul numeric, grafică, aplicații, Editura Teora, București, 2001.

• D. Trif — Tehnici de simulare numerică cu MATLAB, Editura InfoData, Cluj-Napoca, 2007.

• D. Higham, N. Higham — MATLAB Guide, SIAM, Philadelphia, 2000.

4

EVALUARE Condiţii Participare la măcar 14 ore de curs și 14 ore de seminar.

Criterii

Prezentarea și elaborarea unei lucrări practice în

săptămâna a 8-a a semestrului. Cei care, din diferite

motive, nu prezintă lucrarea la timp, vor prezenta

lucrările în săptămâna a 12-a și după aceea se vor putea

prezenta la examen.

Forme

Examen scris și oral, constând din: exerciţii, probleme și

o tematică teoretică, selectată din materia predată la curs

şi comunicată studenţilor în săptămâna a 12-a a

semestrului.

Formula

notei finale

50% va rezulta din lucrări de control, prezenţă şi

activitate la seminarii, 50% va rezulta ca medie a notei la

proba scrisă cu nota la proba orală.

5

1. NOȚIUNI DE BAZĂ

Curs 1

6

1.1. Reprezentarea numerelor reale

Deoarece mașinile de calcul au resurse finite, ele lucrează doar cu o submulțime finită de nume-re reale.

Notă. Datorită modului în care sunt reprezentate, elementele lui se numesc numere în virgulă mobilă.

Orice alt număr real e aproximat prin cel mai apropiat număr în virgulă mobilă: , altfel spus e trunchiat și rotunjit:

7

Cum trunchiază și rotunjește calculatorul numerele

• Să considerăm numărul rațional și reprezentarea sa zecimală infinită periodică: Când cerem programului MATLAB să afișeze acest număr, el o face cu patru zecimale: primele trei sunt exacte, iar a patra e rotunjită prin adaos, din cauză că a cincea este

8


• Să considerăm numărul rațional și reprezentarea sa zecimală infinită periodică: MATLAB afișează acest număr cu pa-tru zecimale exacte; a patra nu mai e rotunjită prin adaos, din cauză că a cincea este

9


• MATLAB — și orice program — trunchiază și rotunjește numerele, dar nu chiar la patru cifre! Ceea ce am văzut e doar un format de afișare (anume short). Reprezentarea internă e mai precisă: cu 16 cifre zecimale. Există și alte formate de afișare a numărului, mai precise decât short.

10

Alte formate de afișare a numerelor în MATLAB

• Comanda ce trebuie utilizată pentru a schimba formatul de afișare este format urmată de numele formatului de afișare dorit:

11

Alte formate de afișare a numerelor în MATLAB

12

Reprezentarea numerelor în calculator (cazul programului MATLAB)

• În general, calculatorul memorează un număr real în modul următor:

• unde este 0 sau 1, este baza de reprezentare (un număr natural mai mare sau egal cu 2), este un număr natural numit mantisă, a cărui lungime t e numărul maxim de cifre memorate , cu Dacă prima cifră neglijată este atunci se adaugă 1 la mantisă. este un număr întreg, numit exponent;

13


• Unicitatea reprezentării

• Mulțimea

• În MATLAB

• Formatul long e cel mai apropiat de reprezentarea internă

• Eroarea relativă e mică:

14


• reprezintă distanța dintre 1 și primul număr în virgulă mobilă mai mare decât el

• În MATLAB și poate fi obținută cu comanda eps.

• Eroarea absolută e mică pentru numerele mici și mare pentru cele mari.

15


• e o mulțime finită, ce nu include numere arbitrar de mari sau arbitrar de mici.

• Să se arate că sunt cel mai mic și cel mai mare număr pozitiv în virgulă mobilă.

• 0 nu e un număr în virgulă mobilă, ci e gestionat separat, la fel ca și Inf sau NaN.

• Numerele în virgulă mobilă sunt mai dense în preajma originii și mai rare în preajma lui ±

16


• În MATLAB

• cele două valori putând fi obținute cu comenzile realmin și realmax

• Distanța între aceste numere și numerele în virgulă mobilă cele mai apropiate de ele este de și respectiv

• Aceste valori ne dau o idee despre densitatea numerelor în virgulă mobilă la cele două extreme ale domeniului reprezentat.

17


• Un număr mai mic decât produce un mesaj de avertisment (underflow) și e tratat apoi ca 0 sau într-un mod special.

• Un număr mai mare decât produce un mesaj de avertisment (overflow) și e memorat apoi în variabila Inf.

• Când calculele conduc la nedeterminări precum 0/0 etc., ele produc ceea ce se numește NaN (not a number).

18

Operații cu numere în virgulă mobilă

• Operațiile algebrice elementare cu numere în virgulă mobilă nu se mai bucură de toate proprietățile operațiilor analoge cu numere reale, deoarece rezultatele sunt trunchiate astfel încât să rămână în domeniul

• Comutativitatea adunării și înmulțirii se păstrează, dar nu și asociativitatea sau distributivitatea!

• Exercițiu. Executați instrucțiunile și interpretați rezultatul! Cu ce valoare a lui b se termină iterația?

19


• În MATLAB comanda eps(a) ne oferă cel mai mic număr care, adunat cu a, dă un rezultat mai mare decât a. Acest număr e cu atât mai mic, cu cât a este mai mic.

• De exemplu, asociativitatea dispare când, grupând termenii într-un anumit fel, se produc erorile de tip overflow sau underflow, dar grupându-i altfel, aceste erori nu se mai produc.

• Exemplu:

20


• În general, se produc erori semnificative de calcul când adunăm numere de semn opus, dar având aproximativ același modul. Spunem că se pierd astfel cifre semnificative.

• Scrieți un program MATLAB care să calculeze pentru o valoare mică a lui x. Interpretați rezultatul.

21


• Scrieți un program MATLAB care să calculeze valoarea funcției

în 401 puncte echidistante din intervalul

• Ind. Puteți folosi linspace

• Reprezentați grafic rezultatul. Care e semnificația lui? Ind. Puteți folosi plot

22

1.2. Erori. Tipuri. Cum le controlăm

• A face erori (a greși) nu e doar omenește! În calculul numeric prezența erorilor este inevitabilă. Ca atare, scopul nostru nu trebuie să fie acela de a le elimina, ci de a le controla încât efectul lor să rămână limitat.

• Erorile apar la mai multe niveluri ale procesului de modelare matematică și numerică a unei realități fizice.

24

Tipuri de erori. Eroarea de modelare

• Astfel, apare mai întâi o eroare de modelare, dată de discrepanța dintre realitatea fizică și modelul matematic — totdeauna imperfect — prin care încercăm să o descriem.

• Ea se traduce în diferența dintre soluția a unei probleme fizice concrete PP și soluția a problemei matematice MP care o modelează.

25

Tipuri de erori. Eroarea de modelare

• Existența unor erori de modelare va restrânge aria de aplicare a modelului matematic la feno-mene de o anumită scară.

– Astfel, mecanica newtoniană e aplicabilă la realitățile vieții cotidiene, teoria relativității e aplicabilă la scară cosmică, iar mecanica cuantică la scară subatomică.

• Eroarea e mai greu de cuantificat. Evaluarea ei depășește cadrul și posibilitățile calculului științi-fic, fiind de competența științelor experimentale.

26

Tipuri de erori. Eroarea de trunchiere

• Cu excepția unor cazuri singulare, soluția a problemei matematice MP nu poate fi găsită în formă explicită.

• Problema matematică e atunci înlocuită cu o problemă numerică NP, a cărei rezolvare furni-zează o soluție aproximativă , dar nu mai presupune decât un șir finit de operații.

• Soluția numerică diferă de soluția teoretică printr-o eroare de trunchiere .

27

Tipuri de erori. Eroarea algoritmică

• Problema e apoi încredințată spre rezolvare calculatorului, ceea ce presupune introduce-rea și propagarea cel puțin a unor erori de rotunjire. Calculatorul furnizează o soluție .

• Ansamblul erorilor care apar prin rezolvarea algoritmică a problemei pe calculator se va numi eroare algoritmică și se va nota cu .

28

Tipuri de erori. Eroarea computațională

• Suma dintre eroarea algoritmică și cea de trunchiere constituie eroarea computațională , cantitatea pe care dorim să o menținem sub control.

• Mai precis, eroarea computațională absolută e

iar eroarea computațională relativă este:

unde este modulul, norma sau altă măsură a mărimii, în funcție de natura lui .

29

Tipuri de erori

30

Controlul erorii computaționale

• Metoda numerică e o aproximare a modelului matematic care se obține, în general, ca o funcție de un parametru de discretizare . – De exemplu, dacă aproximăm integrala cu sumele

Darboux, parametrul de discretizare e norma diviziunii.

– Dacă aproximăm derivata a unei funcții f prin . , parametrul de discretizare este h.

• Metoda numerică se numește convergentă dacă soluția numerică tinde la soluția matematică exactă atunci când .

31


• Metoda numerică se numește convergentă de ordin p dacă, în plus, eroarea computațională (absolută sau relativă), ca funcție de h, e mărginită astfel:

unde C e o constantă, iar p e strict pozitiv.

• Putem, atunci, limita eroarea computațională, alegând h suficient de mic.

32


• Dacă inegalitatea precedentă se poate înlocui cu egalitate (aproximativă), putem pune în evidență ordinul de convergență reprezentând grafic eroarea computațională ca funcție de , pe o scală logaritmică, altfel spus, reprezentând ca funcție de .

• În acest caz, dependența e liniară, iar ordinul de convergență este dat de panta dreptei. Într-adevăr:

• În MATLAB, dacă vectorii x și y conțin abscisele și ordonatele punctelor graficului unei funcții, reprezentarea graficului în scală logaritmică se realizează cu comanda loglog(x,y)

33

Comportamentul erorii, în scală logaritmică, pentru două metode numerice diferite: linia continuă reprezintă o aproximare de ordinul I, cea întreruptă — o aproximare de ordinul al doilea.

34

c

Notă. Eroarea computațională nu e o mărime calculabilă, deoarece depinde de soluția matematică exactă a problemei, care e necunoscută. Din acest motiv se introduc așa numiții estimatori de eroare, niște cantități calculabile care pot fi utilizate pentru a estima eroarea. Estimatorii de eroare nu vor fi discutați deocamdată.

35


1.3. Costuri computaționale

• Costul computațional al unui algoritm este nu-mărul de operații în virgulă mobilă, necesare pentru execuția sa.

• Esențială nu e cunoașterea exactă a numărului de operații executate de un algoritm, ci doar ordinul de mărime al numărului de operații, ca funcție de un parametru d legat de dimensiu-nea problemei.

36

Costuri computaționale

Un algoritm are

• complexitate constantă, dacă cere un număr de operații independent de d, adică operații;

• complexitate liniară, dacă execută operații

• și, mai general, complexitate polinomială, dacă execută operații;

• complexitate exponențială, dacă necesită operații

• și complexitate factorială, dacă necesită operații 37


Exemple • Pentru a calcula prin metoda uzuală produsul între o

matrice pătratică de dimensiune n și o coloană de aceeași dimensiune, sunt necesare n(2n-1) operații. Algoritmul are deci complexitate pătratică.

• Calculul determinantului prin dezvoltarea sa după o linie sau o coloană are complexitate factorială! Unui calculator capabil să execute 1 Peta-flops i-ar trebui 20 de ani să ducă la bun sfârșit calculul pentru o matrice de dimensiune 24.

• Din fericire, există și algoritmi mai eficienți: algoritmul lui Strassen, de exemplu, are o complexitate de

38


Un alt factor relevant pentru eficacitatea unui algo-ritm este timpul necesar de accesare a memoriei calculatorului (care depinde de modul în care algoritmul a fost scris). Acest timp depinde și de parametrii mașinii. Cele mai rapide mașini de astăzi pot executa peste 40 Tera-operații în virgulă mobilă pe secundă (Tera-flops). (1 Tera = 1012) Timpul CPU alocat MATLAB de la deschiderea sa este returnat de comanda cputime, iar timpul scurs între două momente t1 și t2 poate fi calculat în MATLAB cu comanda etime. Astfel, cu ajutorul acestor comenzi, se poate calcula timpul CPU necesar execuției unei secvențe de instrucțiuni sau, în combinație cu clock, care returnează data și ora curente, timpul total de execuție al unui algoritm.

39

1.4. Reprezentarea numerelor complexe în MATLAB

40

2. Ecuații neliniare

Cursul al treilea

41

Ecuații neliniare. Privire generală

Una din problemele frecvent întâlnite în calculul științific este aceea a rezolvării unor ecuații de forma

altfel spus, aceea a determinării zerourilor sau rădăcinilor unei funcții reale .

Pentru funcții liniare, soluția se determină imediat. Dacă însă e o funcție neliniară, soluția problemei nu se poate afla într-un număr finit de pași. Nici măcar pentru polinoame de grad strict mai mare ca 4 nu există formule explicite pentru rădăcini.

42

Ecuații neliniare. Privire generală

Pentru a rezolva problema, se folosesc metode iterative. Pornind de la una sau mai multe valori inițiale, se construiesc succesiv valorile , care converg către o rădăcină a funcției .

Altfel spus, valoarea inițială e o estimare a priori (o ghicire) a rădăcinii (sau valorile inițiale . și sunt limitele unui interval ce conține rădăcina), iar . sunt aproximările din ce în ce mai bune ale rădăcinii, pe care metoda iterativă le produce.

43

Ecuații neliniare. Probleme reprezentative

Problema fondului de investiții. La începutul fiecărui an, clientul unei bănci depozitează euro într-un fond de investiții. După n ani, retrage un capital de M euro. Care a fost dobânda medie r anuală?

Relația dintre parametrii problemei fiind

deducem că r e soluția ecuației

unde 44


Ecuația de stare a unui gaz. Vrem să de-terminăm volumul V pe care îl ocupă un gaz având temperatura T și presiunea p. Pentru aceasta, va trebui să rezolvăm, în V, ecuația (neliniară) de stare a gazului:

În această relație, N este numărul de molecule conținute în volumul V, k este constanta lui Boltzmann, iar a și b sunt doi parametri ce depind de gazul considerat.

45

Metoda înjumătățirii intervalului

• O funcție f, continuă pe intervalul [a, b] și care satisface relația f(a) f(b) < 0, se anulează măcar o dată în intervalul (a, b) — teorema valorilor intermediare.

• Pentru simplitate, vom presupune că f are o sin-gură rădăcină în intervalul (a, b). Dacă există mai multe, putem încerca să reducem intervalul prin examinarea graficului, comanda MATLAB ce poa-te fi folosită fiind fplot. Dar e, probabil, suficient să remarcăm că, după un număr finit de înjumă-tățiri ale intervalului, în intervalul construit rămâ-ne o singură rădăcină (dacă sunt în număr finit).

46


Câteva iterații cu metoda înjumătățirii intervalu-lui pentru funcția .

47


• Pornim deci cu o funcție f și cu un interval [a, b] cu proprietatea că f(a) f(b) < 0 și căutăm o rădăcină în intervalul (a, b).

• Strategia metodei constă în a înjumătăți intervalul dat și în a păstra subintervalul între capetele căru-ia funcția f își schimbă semnul, deoarece acesta va conține cu siguranță o rădăcină. Procedeul se repetă până când lungimea subintervalului devine mai mică decât precizia dorită de aproximare.

48

Metoda înjumătățirii intervalului Dacă notăm cu subintervalul închis ales în etapa , cu mijlocul său și cu un zero al funcției f situat în , avem:

–

– și de unde și deci

Conform teoremei lui Cantor, avem atunci

iar conform lemei cleștelui,

Cum însă iar f e continuă, avem 49


• Astfel constituie aproximări din ce în ce mai bune ale rădăcinii a lui f. Mai precis, eroarea pe care o facem luând pe în loc de este

50

2. Ecuații neliniare (continuare)

Cursul al patrulea

Metode de punct fix

51


Probleme de dinamică a populației. Pentru a studia evoluția în timp a unei populații (de bacterii, de exemplu), considerăm ecuația care stabilește o relație între numărul de indivizi . dintr-o generație și numărul de indivizi din generația următoare. Există mai multe modele teoretice pentru rata de creștere R(x) a populației, cele mai cunoscute fiind:

1. modelul malthusian 2. creșterea cu resurse limitate 3. modelul pradă-prădător

52

1. Modelul lui Thomas Malthus (1766–1834)

2. Modelul creșterii cu resurse limitate (Pierre François Verhulst, 1804–1849) ameliorează modelul lui Malthus considerând că creșterea unei populații e limitată de resursele disponibile:

3. Modelul pradă-prădător e o generalizare a modelului lui Verhulst, care ia în calcul prezența unei populații antagoniste:


53


Dinamica populației este deci complet descrisă de procesul iterativ

unde x(k) reprezintă numărul de indivizi prezenți în generația k , fiind generația ancestrală.

Efectivele staționare (de echilibru) x∗ ale populației considerate vor fi soluțiile ecuației:

altfel spus,

Acesta este un exemplu de problemă de punct fix.

54

Ecuații neliniare. Probleme de punct fix

Definiție. Se numește punct fix al unei funcții . un număr real încât Metoda iterativă de determinare a punctului fix. Când funcția admite un punct fix, acesta poate fi uneori determinat, cu aproximație oricât de bună, cu ajutorul algoritmului iterativ: (1) Într-adevăr, în anumite ipoteze de regularitate, unde oricare ar fi data inițială (sau doar când e aleasă într-o vecinătate a punctului fix).

55


Reprezentarea câtorva iterații de punct fix pentru două funcții de iterație diferite. În stânga, metoda converge (funcția e cos x), în dreapta nu.

56

Experiment

Aplicăm repetat cos valorii reale 1. Obținem succesiv valorile

care tind către

57


Reducerea unor probleme la probleme de punct fix. Metoda lui Newton poate fi privită și ca un algoritm de punct fix, pentru funcția de iterație: În general, iterațiile de punct fix pot fi folosite pentru calculul zerourilor unei funcții, iar funcția fn nu e singura alegere posibilă. O alegere evidentă, dar nu foarte utilă în practică, e Un alt exemplu, mai util, este În fine, un alt exemplu e dat de dacă Metoda înjumătățirii intervalului nu poate fi văzută ca o metodă de punct fix pentru că ea nu ia în calcul valorile funcției.

58

Ecuații neliniare. Probleme de punct fix. Teoreme de convergență

Așa cum am văzut în figură, iterațiile de punct fix pot să nu conveargă. Totuși, 1. Remarcă. Dacă șirul definit în (1) este conver-

gent la , iar e funcție continuă, atunci e un punct fix pentru funcția .

2. Teoremă. O funcție continuă f:[a,b][a,b] posedă măcar un punct fix.

3. Teorema de punct fix a lui Banach. Fie I un sub-interval închis al lui R și f:II o funcție cu proprietatea Atunci f admite un punct fix unic , iar șirul (1) tinde către a, oricum l-am alege pe x0. 59


În practică e adesea dificil să găsim un interval pe care să fie îndeplinite ipotezele din teore-mele precedente. În aceste cazuri va fi util rezulta-tul următor, de convergență locală.

1. Teorema lui Ostrowski. Fie un punct fix al unei funcții de clasă pe o vecinătate convenabilă a lui . Dacă atunci există . încât pentru orice cu șirul definit în (1) converge la . Mai mult, are loc

60


În condițiile ultimelor două teoreme avem asigurată convergența (cel puțin) liniară. Dar metoda lui Newton nu e singura procedură iterativă care se bucură de convergență pătratică. În acest sens avem următoarea completare a teoremei precedente: 1. Teorema de convergență pătratică. În ipotezele

teoremei precedente, presupunând că este de două ori derivabilă în , cu avem

61


În continuare ne propunem să dăm demonstra-ții complete pentru teoremele precedente. Dacă , avem, folosind con-tinuitatea lui , de unde Remarca 1. Pentru Teorema 2, considerăm funcția care verifică și , căci . Cum e continuă, ea are proprietatea lui Darboux, deci se anulează într-un c aflat între a și b:

62

Ecuații neliniare. Probleme de punct

fix. Teoreme de convergență

63



64



65



66



67



68



69


Să discutăm un pic caracterul liniar al convergenței asigurat de teoremele lui Banach și Ostrowski. În cazul teoremei lui Banach, eroarea la pasul k+1 este mai mică decât de L ori eroarea de la pasul k. Putem spune că avem un factor de convergență de cel mult L, cu L < 1. În cazul teoremei lui Ostrowski, pentru k îndeajuns de mare, eroarea la pasul k+1 este de ori eroarea de la pasul k. Putem spune că avem un factor asimptotic de convergență de . Reluând calculele în demonstrația Teoremei lui Ostrowski, vedem că dacă , șirul (1) nu poate converge la punctul fix, iar dacă , nu ne putem pronunța. Putem observa convergență sau divergență, după funcția sau după valoarea inițială aleasă. Exemplu. Să se analizeze funcțiile cos(x) și x2 – 1. 70


• Cu cât factorul de convergență (asimptotic sau maximal) e mai mic, cu atât convergența e mai rapidă.

• În ipotezele ultimei teoreme, putem vorbi de un factor asimptotic de convergență pătratică de : pentru k îndeajuns de mare, eroa-rea la pasul k+1 este de ori pătratul erorii de la pasul k. Și aici, cu cât factorul de convergență e mai mic, cu atât convergența e mai rapidă.

71

Ecuații neliniare. Probleme de punct fix. Oprirea iterațiilor

Ca și la metoda lui Newton, iterațiile vor fi între-rupte când diferența dintre două estimări succesive ale punctului fix nu depășește o limită de toleranță fixată dinainte. Acest fapt poate fi justificat rațional în felul urmă-tor. Din teorema lui Lagrange, găsim punctele între și încât să avem: Dar de unde

72

Ecuații neliniare. Probleme de punct fix. Oprirea iterațiilor

• În consecință, dacă într-o vecinătate a lui , diferența dintre două estimări succesive este un estimator de eroare satisfăcător. Acest fapt are loc când convergența e pătrati-că, de pildă în cazul metodei lui Newton.

• Estimatorul devine nesatisfăcător când deriva-ta se poate apropia oricât de mult de 1.

73

Ecuații neliniare. Algoritm de punct fix. Algoritmul de punct fix e implementat în programul pctfix.m de mai jos. Variabilele de intrare sunt: funcția de iterație fi, estimarea inițială x0, toleranța admisă tol și numărul maxim de iterații nmax. Variabilele de ieșire sunt: estimarea găsită x a punctului fix și numărul necesar de iterații niter. function [x, niter]=pctfix( fi, x0, tol, nmax, varargin) % Calculeaza un punct fix al functiei fi, % pornind de la estimarea initiala x0 if nargin == 2 % valori implicite tol = 1.e-4; nmax = 100; elseif nargin == 3 % valoare implicita nmax = 100; end 74

Ecuații neliniare. Algoritm de punct fix.

x = x0; diff = tol + 1; niter = 0; while niter <= nmax & diff >= tol % Iteratia de punct fix xnou = feval( fi, x, varargin{:}); % x = fi (x) diff = abs(x - xnou); x = xnou; niter = niter + 1; end if niter >= nmax fprintf(’Nu converge in numarul maxim de iteratii admis\n’); end return

75

Ecuații neliniare. Exemplu

Aplicăm algoritmul de punct fix problemei de dinamică a populației enunțate la început, cu valorile parametrilor r = 3 și K = 1. Funcțiile de iterație, pentru modelul lui Verhulst și pentru modelul pradă-prădător, sunt respectiv: Pentru a alege o estimare inițială convenabilă, se pot desena graficele funcțiilor de iterație. Intersecția lor cu prima bisectoare va fi punctul fix căutat!

76

Ecuații neliniare. Exemplu

78

2. Ecuații neliniare (continuare)

Cursul al cincilea

Metoda de accelerare a lui Aitken

Determinarea rădăcinilor polinoamelor

79

Metoda de accelerare a lui Aitken Vom prezenta o tehnică ce permite accelerarea conver-genței unui șir construit printr-o metodă de punct fix.

80


81


Nu vom prezenta demonstrații complete, ci vom da o justificare euristică metodei lui Aitken. În acest scop, să observăm că dacă șirul con-verge liniar la punctul fix, atunci, pentru k mare,

de unde sau

Aceasta ne sugerează să introducem un nou șir de iterații: 82


sau

Această expresie e cunoscută ca formula de extrapolare a lui Aitken. Ea este o nouă metodă de punct fix pentru funcția de iterație:

83


De fapt funcția nu e definită pentru , atât numitorul, cât și numărătorul anulându-se în acest caz. Presupunând derivabilă în cu și aplicând regula lui l’Hôpital, obținem totuși

Putem deci prelungi pe prin continuitate prin 84


• Se poate vedea că, dacă , atunci chiar și atunci când . Cazul are loc atunci când e rădăci-nă multiplă a lui

• Se poate verifica, în general, că punctele fixe ale lui sunt tocmai punctele fixe ale lui .

• Se poate arăta că metoda lui Aitken converge chiar dacă metoda obișnuită diverge, cu condiția să avem .

85

Metoda de accelerare a lui Aitken • Avantajele metodei lui Aitken se pot constata din următoarea

86


Programul aitken.m implementează metoda de accelerare a lui Aitken. Variabilele de intrare sunt: • phi: o funcție introdusă printr-un fișier .m sau ca

funcție inline; • x0: data inițială, de la care începem iterațiile; • tol: toleranța în ce privește distanța între două

estimări succesive; valoare implicită: 1.e-04; • nmax: numărul maxim de iterații permis; valoare

implicită 100; • varargin: lista parametrilor funcției phi. Variabilele returnate sunt x: punctul fix și numărul niter al iterațiilor efectuate.

87


Exemplu. Ne pro-punem să găsim, printr-o metodă de punct fix, rădăcina unică a funcției: .

Vom folosi urmă-toarele funcții de iterație:

90


Cum vedem, con-vergența e extrem de rapidă. Pentru comparație, meto-da punctului fix ce-re cu a doua func-ție 18 iterații și nici măcar nu converge cu prima funcție. Într-adevăr, .

91

Metode de punct fix. Recapitulare 1. Un număr satisfăcând e numit punct

fix al funcției . 2. Pentru calculul său repetăm calculul

până când două aproximații succesive nu mai diferă prin mai mult decât o toleranță prag fixată dinainte.

3. Algoritmul converge când sunt îndeplinite anumite condiții asupra funcției și a derivatei sale. Convergen-ța e îndeobște liniară, dar, când , ea e pă-tratică.

4. Pentru orice funcție de iterație e posibil să construim un nou șir de aproximații succesive cu metoda lui Aitken, care converge mai repede.

5. Algoritmii de punct fix pot fi folosiți și pentru calculul zerourilor unei funcții. 92


Considerăm acum cazul când f este o funcție po-linomială de grad n, altfel spus are forma: Următoarele fapte sunt adevărate despre poli-noame: • Dacă e o rădăcină complexă a unui poli-

nom cu coeficienți reali, atunci și conjugata sa e și ea o rădăcină a polinomului.

• Nu există o formulă explicită de calcul al rădăcini-lor unui polinom generic de grad (Teorema lui Abel).

93

Determinarea rădăcinilor polinoamelor Fapte:

• Fie k numărul rădăcinilor reale pozitive, fiecare socotită cu multiplicitatea ei, ale unui polinom de grad n. Fie s numărul schimbărilor de semn ale coeficienților săi. Atunci și este par (Regula semnelor a lui Descartes).

• Toate rădăcinile unui polinom sunt conținute în discul din centrat în origine, de rază , unde (Teorema lui Cauchy).

Aceste ultime două reguli ne pot ajuta să alegem o estimare inițială a rădăcinilor căutate.

94


Exemplu. Polinomul de gradul 6

are rădăcinile , deci are o singură ră-dăcină reală pozitivă; el are 5 schimbări de semn. Se îndeplinesc concluziile din Teorema lui Descartes. Teorema lui Cauchy ne asigură că toa-te rădăcinile sunt de normă cel mult egală cu 9. Evident, suntem interesați de estimări mai bune decât cea dată de Teorema lui Cauchy.

95


Algoritmul lui Hörner permite evaluarea poli-nomului (și a derivatei sale) într-un punct dat, precum și determinarea câtului și restului la împărțirea polinomului prin

Metoda reducerii gradului. Vom vedea apoi că acest algoritm se poate combina convenabil cu metoda lui Newton de găsire a unui rădăcini, pentru a crea un algoritm de determinare a tuturor rădăcinilor unui polinom, numit metoda reducerii gradului.

96

Determinarea rădăcinilor polinoamelor Putem scrie o funcție polinomială în forma:

Să observăm că pentru a evalua această expresie sunt necesare n adunări și n înmulțiri, în vreme ce, în forma inițială, ar fi necesare n adunări și 2n-1 înmulțiri (de ce doar 2n-1?). Acesta e numărul de operații efectuate de algoritmul împărțirii sintetice:

97


în care, evident, Considerând polino-mul (pe care-l vom numi asociat lui f)

acesta are proprietatea remarcabilă de a fi toc-mai câtul împărțirii lui f(x) la x-z, în vreme ce este chiar restul aceleiași împărțiri. Altfel spus,

Verificați acest lucru prin identificarea coeficien-ților! z e rădăcină dacă și numai dacă

98


Să implementăm algoritmul împărțirii sintetice într-un program horner.m.

Variabilele de intrare vor fi a și z. Vectorul n-di-mensional a memorează coeficienții lui f, în ordine inversată, iar z memorează argumentul în care dorim să evaluăm polinomul.

Variabilele de ieșire vor fi b, care memorează, în ordine inversată, coeficienții câtului și r, restul la împărțirea cu .

99

Determinarea rădăcinilor polinoamelor. Metoda Newton-Hörner

Pentru a calcula toate rădăcinile unui polinom vom adopta următoarea strategie de reducere a gradului: Cât timp gradul polinomului nu e 0,

1. găsim o rădăcină z printr-o metodă adecvată de aproximare; vom folosi metoda lui Newton;

2. reducem gradul polinomului, împărțindu-l la x – z. Folosim algoritmul lui Hörner.

101


• Avantajul acestei metode constă în aceea că metoda lui Newton se cuplează convenabil cu schema lui Hörner.

• Într-adevăr, notând polinomul la etapa n (când are gradul n) cu , avem

de unde

102


• Metoda lui Newton capătă forma

în care și , și coeficienții lui , și valoarea lui sunt evaluați cu schema lui Hörner.

• După atingerea preciziei dorite în determina-rea rădăcinii, reducem gradul tot cu schema lui Hörner.

103


• Continuăm cu aproximarea unei rădăcini a po-linomului de grad mai mic și așa mai departe, până când toate cele n rădăcini au fost calcula-te.

• Observăm că, pentru a obține și rădăcinile complexe r ∈ C \ R, trebuie aleasă o estimare inițială cu partea imaginară nenulă.

104


• Metoda Newton-Hörner e implementată în programul următor. Folosim, în mod implicit, un număr maxim de 100 de iterații și tolerăm o discrepanță maximă de 1.e-04 între două aproximații succesive. Ca și în programul horner.m, coeficienții polinomului inițial sunt stocați în vectorul a, în ordine inversată.

• Programul returnează în vectorii roots și iter rădăcinile calculate și numărul de iterații care au fost necesare pentru a atinge precizia dată de tol.

105


• Observație. Pentru a minimiza propagarea erorilor de rotunjire, în procesul de reducere a gradului este recomandabil să aproximăm mai întâi rădăcina r1 cu modulul cel mai mic, și abia apoi să trecem la aproximarea celorlalte rădăcini r2 , r3 , ... până se ajunge la cea cu modulul cel mai mare.

• Din acest motiv, estimarea inițială se ia adesea r0=0.

108


Exemplul 1. Calculul rădăcinilor polinomului

Rădăcinile sale sunt

Metoda converge repede, iar rezultatele obținute sunt precise.

109


Exemplul 2. Calculul rădăcinilor polinomului

Acesta are pe 1 drept rădăcină triplă și pe 4 drept rădăcină simplă.

E una din situațiile când erorile de rotunjire se propagă. Metoda nu mai e la fel de eficientă. 110

COMENZI MATLAB SPECIALIZATE

111

Comenzi Matlab specializate Cele mai sofisticate metode de calcul al zero-urilor unei funcții neliniare sunt combinații ale diverșilor algoritmi, în încercarea de a obține o convergență mai bună către soluție.

De exemplu, funcția Matlab predefinită fzero calcu-lează zerourile unei funcții folosind așa-numita metodă Dekker-Brent. Această comandă are sintaxa de bază

fzero (fun, x0)

unde argumentul fun e o funcție inline sau definită într-un fișier .m, iar x0 este o estimare inițială.

112

Comenzi Matlab specializate

Să rezolvăm problema fondului de investiții: function y= fond(r)

y=6000 - 1000*(1+r)/r*((1+r)ˆ5 - 1);

=================================== ≫ [I, res, flag] = fzero (’fond’, x0) Zero found in the interval: [0.028471, 0.492]. I = 0.06140241153653 res = -1.818989403545857e-012 flag = 1 Dacă variabila de ieșire flag este < 0, înseamnă că funcția fzero nu poate găsi rădăcina.

113

Comenzi Matlab specializate Pentru calculul rădăcinilor unui polinom, în afară de metoda Newton–Hörner, există numeroase alte metode de aproximare. Astfel, comanda Matlab roots folosește o metodă prin care se caracterizează zerourile unei funcții ca și valorile proprii ale unei matrice speciale asociate. Pentru rezolvarea sistemelor de ecuații neliniare, există generalizări ale metodei lui Newton, numite metode quasi–Newton. Comanda Matlab zero = fsolve(’fun’, x0) permite aflarea unei soluții pentru sisteme de ecuații. Variabi-la fun este numele unui fișier .m în care a fost definită forma funcțiilor iar x0 este estimarea inițială a soluției.

114

Comenzi Matlab specializate De exemplu, pentru rezolvarea sistemului : scriem instrucțiunile: function fx = expsist(x) fx(1) = x(1)�2+x(2)�2-1; fx(2) = sin(pi*0.5*x(1))+x(2)�3; ==================================== » x0 = [1, 1]; » sol = fsolve(’expsist’, x0) sol = 0.4761 -0.8794 Pentru a obține soluția cu semn schimbat, pornim de la estimarea inițială x0=[-1,-1].

115

Curs 6

Aproximarea funcțiilor și a datelor. Interpolare

116

Argument A aproxima o funcție f înseamnă a o înlocui cu o funcție mai simplă 𝑓 și care i se poate substitui acceptabil în calcule. Această strategie este folosită de exemplu în integrarea numerică, unde înlocuim funcția de integrat cu o funcție a cărei integrală e mai ușor de calculat (de exemplu cu o funcție polinomială). Vom aborda această problemă într-un curs ulterior. În alte împrejurări, de exemplu atunci când provine din măsurători, funcția 𝑓 poate să nu fie disponibilă decât parțial, prin valorile (aproximative) pe care le ia în câteva puncte, aflate în număr finit. Dorim atunci să construim o funcție care ar putea reprezenta legea empirică aflată în spatele datelor experimentale (și care le-ar putea rezuma într-o formulă simplă). Sa considerăm câteva exemple tipice care ilustrează această abordare. 117

Probleme tipice

Problema 1 (Climatologie) — Temperatura aerului în apropierea solului depinde de concentrația 𝐾 de acid car-bonic (H2CO3) din aer. În tabelul următor (luat din Philo-sophical Magazine 41, 237 (1896)) e dată diferența între temperatura medie observată pentru diferite latitu-dini și pentru patru valori diferite ale concentrației 𝐾 și temperatura medie corespunzătoare unei valori de refe-rință 𝐾 . Drept concentrație de referință 𝐾 e luată valoa-rea măsurată în 1896, normalizată la 1. Se pune problema găsirii unei funcții care, pe baza datelor existente, să furnizeze o valoare aproximativă a temperaturii medii la orice latitudine și la orice concen-trație 𝐾 de acid carbonic.

118

Probleme tipice. Climatologie

119

Probleme tipice

Problema 2 (Finanțe) — În figura următoare trasăm evo-luția prețului unei acțiuni la bursa din Zürich vreme de doi ani. Curba a fost obținută unind cu linii drepte cotațiile zilnice la ora închiderii. Această reprezentare presupune implicit că prețurile evoluează liniar în cursul unei zile. (Această aproximare se va numi interpolare afină compo-zită). Se pune întrebarea dacă se poate prezice, folosind acest grafic, evoluția prețului acțiunii pe un scurt interval de timp după ultima cotație înregistrată. Problema se va rezolva cu o tehnică de aproximare numită metoda celor mai mici pătrate.

120

Probleme tipice. Finanțe

121

Probleme tipice

Problema 3 (Biomecanică) — Considerăm un test mecanic realizat pentru a stabili legătura dintre efortul 𝜎 la care e supus (măsurat în Megapascali; 1 MPa = 100 N/cm2) și deformația 𝜖 pe care o suferă un eșantion de țesut biologic (un disc intervertebral, v. figura).

Pornind de la datele înregistrate în tabel, vom estima deformația corespunzătoare unui efort de σ = 0.9 MPa.

122

Probleme tipice. Biomecanică

123

Aproximarea funcțiilor prin polinoame Taylor

O funcție 𝑓: 𝐼 → derivabilă până la ordinul n într-un punct 𝑥0 ∈ 𝐼 poate fi aproximată prin polinomul Taylor ce-i este asociat:

𝑇𝑛𝑓 𝑥 ≔ 𝑓 𝑥0 +1

1!𝑓′ 𝑥0 𝑥 − 𝑥0 +

1

2!𝑓′′ 𝑥0 𝑥 − 𝑥0

2 + ⋯+1

𝑛!𝑓(𝑛)(𝑥0) 𝑥 − 𝑥0

𝑛

Aproximarea este bună într-o vecinătate a punctului 𝑥0. Mai precis, restul 𝑅𝑛𝑓 ≔ 𝑓 − 𝑇𝑛𝑓 poate fi pus în forma lui Peano:

𝑅𝑛𝑓(𝑥) =𝛼(𝑥)

𝑛!(𝑥 − 𝑥0)𝑛

cu 𝛼: 𝐼 → continuă și nulă în 𝑥0.

124


Mai mult, dacă avem derivabilitate de ordinul 𝑛 + 1 pe un întreg interval, atunci restul poate fi pus în forma lui Lagrange:

Teoremă (Formula lui Taylor cu restul lui Lagrange) — În ipotezele deja asumate pentru 𝑓, presupunând că 𝑥 ∈ 𝐼 și că 𝑓 e derivabilă de 𝑛 ori pe întreg intervalul [𝑥0, 𝑥] (sau [𝑥, 𝑥0]) și de 𝑛 + 1 ori pe intervalul (𝑥0, 𝑥) (sau (𝑥, 𝑥0) ), există

θ ∈ [𝑥0, 𝑥], încât 𝑅𝑛𝑓(𝑥) =𝑓(𝑛+1)(𝜃)

(𝑛+1)!(𝑥 − 𝑥0)𝑛+1.

125


Această tehnică e costisitoare din punct de vedere computațional, pentru că ea cere cunoașterea funcției f și a derivatelor sale până la ordinul n (egal cu gradul polinomului) într-un punct dat 𝑥0. Mai mult, polinomul Taylor s-ar putea să nu mai reprezinte cu acuratețe funcția dacă ne îndepărtăm prea mult de punctul 𝑥0. Urmăriți, pe graficele următoare, îmbunătățirea aproximării și validitatea ei pe un interval mai lung, odată cu creșterea ordinului n, dar și îndepărtarea bruscă a graficului polinomului aproximant de cel al funcției aproximate, de îndată ce am ieșit dintr-o anumită vecinătate a punctului 𝑥0.

126


127


128


Alt exemplu: Considerăm funcția f(x) = 1/x pe intervalul [1,3] și polinomul lui Taylor de grad n = 10 în punctul x0 = 1. După cum se vede din figura următoare, aproximația este foarte bună într-o mică vecinătate a punctului x0 = 1, apoi devine nesatisfăcătoare. (Linia continuă este f(x), iar linia punctată este T10 f(x).) Figura arată și interfața grafică a funcției MATLAB taylortool, prin care se pot calcula și reprezenta grafic polinoamele Taylor de grad N asociate unei funcții date f(x), într-un punct a și pe un interval oarecare. Se apelează prin: ≫ taylortool

129


130

Alte abordări

În continuare vom descrie metode de aproximare bazate pe abordări alternative: interpolare, interpolare compozită, metoda celor mai mici pătrate.

131

Interpolare

132

Curs 7

Derivare și integrare numerică

133

Derivare și integrare numerică În acest curs vom descrie metode de aproximare numerică a derivatei și a integralei unei funcții. Argument. În ce privește integrala, este cunoscut că primitivele unor funcții elementare nu sunt funcții elementare și, deci, nu pot fi găsite în formă explicită. Uneori, chiar și a-tunci când dispunem de o formulă exactă pentru primitivă (sub forma unei serii, de exemplu), evaluarea acesteia e o pro-blemă la fel de dificilă din punct de vedere numeric. Alteori funcția de derivat sau de integrat e definită ca soluție a unei ecuații diferențiale, ce nu poate fi rezolvată explicit, sau nu e cunoscută decât sub formă de tabel, în câteva puncte (de exemplu atunci când apare ca rezultat al unor măsurători). În toate aceste cazuri e necesar să luăm în considerare metode numerice, în vederea obținerii unei aproximări a can-tității care ne interesează, indiferent dacă funcția e greu sau ușor de derivat sau de integrat.

134

Probleme tipice

Problema 1 (Hidraulică) Un lichid se scurge dintr-un vas ci-lindric de rază 𝑅 = 1 𝑚 printr-un orificiu de rază 𝑟 = 0,1 𝑚

135

Se cere determinarea vitezei de golire 𝑞′(𝑡) și compararea ei cu viteza teoretică prezisă de legea lui Torricelli

𝑞′ 𝑡 = −𝛾(𝑟/𝑅)2 2𝑔𝑞(𝑡), unde 𝑔 e accelerația gravitațională, iar 𝛾 = 0,6 e un factor de corecție.

Să prezentăm câteva probleme tipice care presupun deriva-rea și integrarea numerică.

aflat la baza vasului. Înălțimea coloanei de li-chid la momentul 𝑡 se notează cu 𝑞(𝑡) și ea se măsoară la fiecare 5 𝑠,obținându-se rezultate-le din tabelul următor.

Probleme tipice Problema 2 (Demografie) — Considerăm o populație formată dintr-un număr foarte mare 𝑀 de indivizi. Re- partiția 𝑁(ℎ) a înălțimilor acestora poate fi reprezen-tată printr-o funcție clopot, caracterizată de media ℎ și de deviația standard 𝜎: Atunci reprezintă numărul de indivizi cu înălțimea cuprinsă între ℎ și ℎ + ∆ℎ. Un exemplu concret se poate urmări în figura următoare. 136

Pentru 𝑀 = 200, ℎ = 1,7 𝑚 și 𝜎 = 0,1 𝑚, a calcula numărul de oameni cu înălțimea cuprinsă între 1,8 și 1,9 𝑚 revine la a calcula aria porțiunii colorate, sub forma unei integrale. 137

Aproximarea derivatelor prin diferențe finite

Să considerăm o funcție derivabilă într-un punct 𝑥 ∈ (𝑎, 𝑏). Dată fiind definiția de-rivabilității lui 𝑓 în 𝑥 , cantitățile

și

sunt aproximări oricât de bune ale lui 𝑓′(𝑥 ) pen-tru ℎ > 0 suficient de mic. Ele sunt numite dife-rențe finite, progresivă și, respectiv, regresivă.

138


O altă idee e aceea de a considera diferența finită centrată

ca aproximare a derivatei. Într-adevăr

139


De fapt aceste cantități aproxi-mează pe 𝑓′(𝑥 ) (adică panta tan-gentei în (𝑥 , 𝑓 𝑥 )) cu panta

• secantei ce trece prin punctele (𝑥 , 𝑓 𝑥 ) și (𝑥 + ℎ, 𝑓 𝑥 + ℎ ),

• a secantei ce trece prin punctele (𝑥 , 𝑓 𝑥 ) și (𝑥 − ℎ, 𝑓 𝑥 − ℎ ),

• respectiv a coardei ce trece prin punctele (𝑥 + ℎ, 𝑓 𝑥 + ℎ ) și (𝑥 − ℎ, 𝑓 𝑥 − ℎ ).

140


Dacă presupunem în plus că 𝑓 este 𝐶2((𝑎, 𝑏)), atunci putem estima eroarea pe care o facem înlocuind derivata cu diferen-ța finită și viteza de convergență a diferenței finite la derivată când ℎ → 0. Într-adevăr, din formula lui Taylor avem: și, respectiv, de unde obținem și adică diferențele finite progresivă și regresivă sunt aproximări de ordinul întâi în raport cu ℎ ale derivatei. 141


Dacă presupunem acum că 𝑓 este 𝐶3((𝑎, 𝑏)), atunci putem dovedi că diferența finită centrată e o aproximare de ordinul al doilea în raport cu ℎ a derivatei când ℎ → 0. Într-adevăr, din formula lui Taylor avem: și de unde obținem, prin scăderea celor două ecuații și împărțire prin 2ℎ: cu și 142


Implementare concretă. Putem în general presupune că valorile funcției 𝑓 sunt disponibile în 𝑛 + 1 puncte echidistante 𝑥𝑖 = 𝑥0 + 𝑖ℎ, 𝑖 = 0,1, … , 𝑛. Valorile 𝑓′(𝑥𝑖) în aceste puncte ale derivatei se aproximează cu ajutorul diferențelor finite în care apar doar valorile funcției 𝑓 în punctele 𝑥𝑖. Să observăm că diferențele centrate nu au sens în capetele 𝑥0 și 𝑥𝑛. Ele se pot înlocui cu estimările

în 𝑥0

în 𝑥𝑛

care se obțin (ca și diferențele finite, de altfel) prin derivarea polinomului de grad 2 ce interpolează 𝑓 în nodurile 𝑥0, 𝑥1, 𝑥2 (resp. 𝑥𝑛−2, 𝑥𝑛−1, 𝑥𝑛)

143


Soluția problemei de hidraulică. Rezolvăm problema pen-tru un pas ℎ = 5. Vom aproxima funcția 𝑞′ 𝑡 în punctele 𝑡 = 0, 5, 10, 15, 20.

144


Trecem rezultatele în următorul tabel, unde ele pot fi comparate și cu rezultatele teoretice așteptate, date de legea lui Torricelli.

Se observă că există acord între rezultatele numerico-experimentale și cele teoretice, aproximarea prin diferențe centrate având cea mai mare acuratețe.

145

Integrare numerică

În această secțiune vom introduce metode numerice de aproximare a integralei Vom considera mai întâi intervalul de integrare închis și mărginit, iar 𝑓 va fi o funcție continuă pe [𝑎, 𝑏]. Vom începe cu câteva formule simple, care sunt cazuri particulare ale formulelor Newton-Cotes. Apoi vom introduce așa-numitele metode gaussiene, care realizează gradul maxim de exactitate pentru un număr dat de evaluări ale funcției 𝑓.

146

Integrare numerică

Orice formulă explicită 𝐼𝑛(𝑓) care furnizează o a-proximare pentru 𝐼(𝑓) e numită formulă de cuadratură sau formulă de integrare numerică. Formula de cuadratu-ră depinde de un parametru n și, în general, aproximează integrala din ce în ce mai bine pe măsură ce n crește. Ideal, ar trebui să avem 𝐼𝑛 𝑓 → 𝐼 𝑓 (𝑛 → ∞). Chiar și atunci când acest lucru nu va fi imediat vizi-bil, practic toate formulele de cuadratură pe care le vom studia vor fi de forma

(1) 𝐼𝑛 𝑓 = 𝛼𝑖𝑓(𝑥𝑖)𝑛𝑖=0 ,

rezultând deci din 𝑛 + 1 evaluări ale funcției 𝑓. Punctele distincte 𝑥0, 𝑥1, … , 𝑥𝑛 se numesc nodurile cuadraturii, iar scalarii 𝛼0, 𝛼1, … , 𝛼𝑛 se numesc coeficienții (sau ponderi-le) cuadraturii.

147

Criterii de evaluare a aproximării

Problema de bază a integrării numerice va fi atunci alegerea nodurilor și a coeficienților astfel încât 𝐼𝑛(𝑓) să fie o bună aproximare a lui 𝐼(𝑓) pentru o clasă largă de funcții 𝑓.

Diverse criterii pot fi utilizate pentru a evalua calitatea aproximării printr-o cuadratură:

1. Eroarea cuadraturii este 𝐸𝑛 𝑓 = 𝐼 𝑓 − 𝐼𝑛(𝑓).

2. Gradul de exactitate al unei cuadraturi e cel mai mare întreg 𝑟 ≥ 0 pentru care

148

Gradul de exactitate Propoziție. Echivalent, gradul de exactitate este acel 𝑟 ∈ pentru care,

𝐸𝑛 𝑥𝑘 = 0, ∀𝑘 ≤ 𝑟

𝐸𝑛 𝑥𝑟+1 ≠ 0. Demonstrație. Într-adevăr, atât integrala exactă, cât și cuadratu-rile de forma (1), sunt operatori liniari, în sensul că 𝐼 𝛼𝑓 + 𝛽𝑔 =𝛼𝐼 𝑓 + 𝛽𝐼(𝑔) și 𝐼𝑛 𝛼𝑓 + 𝛽𝑔 = 𝛼𝐼𝑛 𝑓 + 𝛽𝐼𝑛(𝑔). Teoremă (Estimarea erorii). Dacă gradul de exactitate al cuadraturii este 𝑟 , atunci pentru orice funcție 𝑓 de clasă 𝐶𝑟+1[𝑎, 𝑏], eroarea cuadraturii este de forma

𝐸𝑛 𝑓 =1

𝑟! 𝑓 𝑟+1 𝑦 𝐸𝑛 𝑥 − 𝑦 +

𝑟 𝑑𝑦𝑏

𝑎

Demonstrație. Se exprimă 𝑓 (în jurul lui 𝑎, de exemplu) cu aju-torul formulei lui Taylor cu rest în formă integrală. Se aplică am-bilor membri ai egalității operatorul 𝐸𝑛 și se ține cont de gradul de exactitate pentru a anula cei mai mulți termeni. În ultimul termen se intervertesc integralele, respectiv suma și integrala. 149

Cuadratură interpolatorie

Practic toate formulele de cuadratură stu-diate de noi se pot obține și înlocuind, în expre-sia integralei, funcția de integrat 𝑓 cu o aproxi-mare a sa 𝑓𝑛 mai ușor de integrat (de exemplu polinomială sau polinomială pe porțiuni) și cal-culând 𝐼(𝑓𝑛) în loc de 𝐼(𝑓). Vom scrie atunci . Vom avea deci:

(2) 150

Cuadratură interpolatorie. Estimarea erorii

Dacă funcția 𝑓 și funcția aproximantă (sau măcar diferența |𝑓 − 𝑓𝑛|) sunt continue pe intervalul [𝑎, 𝑏] , atunci eroarea cuadraturii 𝐸𝑛 𝑓 =𝐼 𝑓 − 𝐼𝑛(𝑓) satisface

Deci, dacă, pentru un anumit 𝑛,

atunci

151

Cuadratură interpolatorie O abordare naturală constă în a alege drept 𝑓𝑛 polinomul de interpolare Lagrange peste un set de 𝑛 + 1 noduri dis-tincte {𝑥𝑖|𝑖 = 0,1, … , 𝑛}. Obținem atunci formula (zisă de cuadratură a lui Lagrange) unde 𝑙𝑖 este polinomul caracteristic Lagrange de grad 𝑛 a-sociat nodului 𝑥𝑖. Gradul de exactitate al cuadraturii lui Lagrange e cel puțin 𝑛. Există mai multe rețete de alegere a nodurilor. Dacă nodurile sunt alese echidistant, obținem așa-numitele formule Newton-Cotes. Dacă nodurile de interpolare sunt alese inteligent, acest grad se poate mări până la 2𝑛 + 1 (în cazul așa-numi-telor formule de cuadratură gaussiene). În continuare, vom considera trei cazuri remarcabile (prin simplitate) de cuadraturi Newton-Cotes, anume pentru 𝑛 = 0, 1 și 2. 152

Integrare numerică Formula dreptunghiurilor (a mijloacelor) e formula de cuadratură a lui Lagrange care se obține luând ca unic nod de interpolare mijlocul intervalului, altfel spus înlocuind funcția originală cu funcția constant egală cu 𝑓

𝑎+𝑏

2. Obținem ca aproximare a

integralei: Dacă 𝑓 e de clasă 𝐶2 pe intervalul (𝑎, 𝑏), eroarea de aproximare e dată de unde e un punct din intervalul (𝑎, 𝑏). Gradul său de exactitate e deci 1. Eroarea poate fi mare.

153

Integrare numerică

Eroarea poate fi mare dacă lungimea intervalului este mare. Acest defect este comun tuturor cuadraturilor de acest tip și el poate fi înlăturat înlocuindu-le cu cuadraturile compozite corespunzătoare. Acestea se obțin împărțind intervalul de integrare într-o diviziune (echidistantă) și aplicând metoda de cuadratură pe fiecare subinterval.

154

Integrare numerică Un procedeu simplu de aproximare a lui I(f) se obține deci partiționând intervalul [a, b] în subintervalele cu și Deoarece putem aplica pe fiecare subinterval 𝐼𝑘 metoda precedentă, aproximând pe 𝑓 cu constanta egală cu valoarea pe care 𝑓 o ia în mijlocul 𝑥 𝑘 al intervalului. Obținem astfel formula compozită de cuadratură (zisă prin dreptunghiuri sau a mijloacelor)

155

Interpretarea geometrică a formulelor de cuadratură prin dreptunghiuri simplă și compozită.

156

Integrare numerică

Dacă 𝑓 e de clasă 𝐶2 pe intervalul [𝑎, 𝑏], eroarea de aproximare e dată de unde e un punct din intervalul [𝑎, 𝑏]. Demonstrația acestei formule se face folosind dezvoltarea Taylor de ordin 2 pe fiecare subinterval și teorema de medie pentru integrale și pentru sume. Cuadratura prin dreptunghiuri are deci o eroare de ordinul doi în raport cu 𝐻. Gradul său de exactitate e 1.

157

Integrare numerică Formula de cuadratură compozită a dreptunghiurilor e implementată în programul Matlab midpointc.m. Parametrii de intrare sunt capetele intervalului 𝑎 şi 𝑏, numărul de subintervale 𝑀 şi funcţia de integrat f.

158

Elemente de Simulare Numerica v3.1

Documents

Transcript of Elemente de Simulare Numerica v3.1