CUPRINS

Introducere 4

Capitolul 1 Erorile de calcul numeric 61.1. Surse de erori 61.2. Propagarea erorilor de calcul 71.3. Algoritmi şi complexitate de calcul 91.4. Metode de programare 11

Capitolul 2 Rezolvarea numerică a ecuaţiilor şi sistemelor de ecuaţii algebrice neliniare 132.1. Introducere 132.2. Rezolvarea numerică a ecuaţiilor neliniare 13

2.2.1. Metoda aproximaţiilor succesive 132.1.1. Metoda Lagrange 152.2.2. Metoda Newton 162.2.3. O teoremă de punct fix 172.2.4. Ordinul metodei 222.2.5. Accelerarea convergenţei 24

2.2.5.1. Metoda Aitken 242.1.1.1. Metoda Steffensen 26

2.2.6. Metoda poziţiei false 282.2.7. Principiul dihotomiei 30

2.3. Rezolvarea numerică a sistemelor neliniare 312.3.1. Metoda aproximaţiilor succesive 322.3.2. Metoda Newton 34

2.4. Exerciţii 35

Capitolul 3 Rezolvarea sistemelor algebrice liniare 363.1. Introducere 363.2. Metode directe 36

3.2.1. Metoda de eliminare a lui Gauss 363.1.1. Factorizarea LU 403.2.2. Factorizarea Cholesky 42

1

3.3. Metode iterative 433.3.1. Metoda iterativă Jacobi 443.3.2. Metoda iterativă Gauss-Seidel 463.3.3. Metoda relaxării 48

3.4. Exerciţii

Capitolul 4 Rezolvarea numerică a problemelor algebrice de valori şi vectori proprii 494.1. Abordarea elementară a subiectului 494.2. Aspecte teoretice generale 50

4.2.1. Clase speciale de matrici 514.1.1. Punerea corectă a problemei 52

4.3. Metoda lui Jacobi 544.4. Probleme de valori proprii generalizate 604.5. Exerciţii 61

Capitolul 5 Aproximarea funcţiilor prin polinoame 625.1. Introducere 625.2. Aproximarea prin interpolare 62

5.2.1. Interpolarea polinomială Lagrange 635.1.1. Algoritmul Aitken 655.2.2. Evaluarea restului la interpolarea Lagrange 675.2.3. Diferenţe divizate 695.2.4. Formula lui Newton de interpolare 705.2.5. Diferenţe finite 725.2.6. Formule de interpolare pe noduri echidistante 725.2.7. Interpolarea polinomială Hermite 745.2.8. Interpolarea prin funcţii spline 76

5.3. Aproximarea în sensul celor mai mici pătrate 785.3.1. Problema fundamentală a aproximării liniare 785.3.2. Teoreme fundamentale 795.3.3. Aproximarea discretă în sensul celor mai mici pătrate 80

5.4. Exerciţii 82

Capitolul 6 Rezolvarea ecuaţiilor diferenţiale ordinare de ordinul I 836.1. Introducere 83

2

6.2. Metode 846.2.1. Metoda Euler 846.2.1. Metoda Euler backward 84

6.3. Generalizări 856.3.1. Metode Runge – Kutta 85

6.3.1.1. Metoda clasică Runge – Kutta de ordin 4 856.3.1.2. Metode Runge – Kutta explicite 866.3.1.3. Metode Runge – Kutta ajustative 876.3.1.4.Metode Runge – Kutta implicite 87

6.3.2. Caracteristici 876.4. Metode alternative 88

Capitolul 7 Integrarea numerică 907.1. Introducere 907.2. Integrarea în puncte echidistante 95

7.2.1. Formulele Newton – Cotes 987.2.1. Metoda dreptunghiurilor 987.2.2. Regula trapezului 997.2.3. Metoda Romberg 1007.2.4. Regula Simpson 1017.2.5. Metoda ajustativă Simpson 103

7.3. Integrarea în puncte neechidistante 1047.3.1. Cvadratura gaussiană 1087.3.1. Cvadratura tanh – sinh 1087.3.2. Cvadratura Clenshaw – Curtis 1097.3.3. Cvadratura Fejer 1127.3.4. Cvadratura ajustativă 112

7.4. Integrarea cu funcţii pondere 1147.5. Metoda Nystrom 1157.6. Formula Euler – MacLaurin 1167.7. T - integrarea 121

Bibliografie 122

INTRODUCERE

3

Ultimele decenii au fost marcate de progresul mijloacelor de calcul. Asistăm la o competiţie între dezvoltarea tehnologică şi dezvoltarea aplicaţiilor, în particular, a celor numerice. Tehnica de calcul a devenit accesibilă pentru categorii tot mai largi de utilizatori. Globalizarea accesului la magistralele informaţiilor organizate în reţeaua Internet a dat o nouă dimensiune utilizării calculatoarelor, revoluţionând domenii întregi de activitate.

Obiectul calculului numeric îl reprezintă găsirea unor metode de aproximare eficientă a soluţiilor problemelor care pot fi exprimate prin modele matematice, eficienţă ce depinde de precizia cerută pentru rezultate şi de uşurinţa implementării. Calculul numeric este una dintre disciplinele matematice ce depinde în cea mai mare măsură de calculatorul numeric.

Drumul parcurs pentru rezolvarea unei probleme dintr-un domeniu oarecare cu ajutorul calculatorului constă în: stabilirea unui model matematic al problemei concrete (model ce se poate încadra într-o categorie cum ar fi: o ecuaţie neliniară, un sistem de ecuaţii liniare sau neliniare), care fiind de multe ori de natură continuă trebuie discretizat; soluţia problemei discretizate trebuie să fie consistentă şi stabilă (robustă); modelul discretizat trebuie transpus într-un algoritm realizabil şi eficient, descris de obicei într-un limbaj de programare evoluat.

Calculul numeric operând cu mărimi variate presupune folosirea tipului real a cărui reprezentare în calculator este aproximativă, apărând erori de rotunjire care se propagă. Deci, o metodă numerică trebuie aleasă ţinând seama de convergenţă, stabilitate, propagarea erorilor şi de analiza complexităţii algoritmului asociat.

Pentru parcurgerea şi utilizarea unui asemenea material, cititorul are nevoie de cunoştinţe de matematică la îndemâna studenţilor care au promovat primul an de studiu al oricărei facultăţi cu profil tehnic, matematico-informatic sau economic.

Metodele numerice sunt prezentate în detaliu, prin discutarea aspectelor de ordin strict matematic şi descrierea algoritmilor cu ajutorul unui limbaj de tip pseudocod.

Lucrarea „Calcul numeric” are şapte capitole.Primul capitolul are un caracter eterogen – la început se prezintă

sursele de erori şi propagarea lor, apoi algoritmi şi complexitate de calcul, iar în final, metode de programare.

Capitolul al doilea are ca obiect rezolvarea numerică a ecuaţiilor şi sistemelor de ecuaţii algebrice neliniare. Sunt prezentate metode de localizare a soluţiei, de aproximaţii succesive şi de accelerare a convergenţei pentru ecuaţii neliniare, precum şi metode numerice de rezolvare a sistemelor algebrice neliniare.

4

Capitolul al treilea este dedicat rezolvării numerice a sistemelor algebrice liniare. Sunt examinate metode directe bazate pe factorizarea gaussiană, precum şi metode de aproximare.

În capitolul al patrulea se prezintă metode de tip Jacobi de rezolvare numerică a problemelor algebrice de valori şi vectori proprii, precum şi generalizarea lor.

În capitolul al cincilea se prezintă aproximarea funcţiilor prin interpolare de tip Lagrange, Hermite şi prin funcţii spline, precum şi aproximarea în sensul celor mai mici pătrate.

În cel de al şaselea capitol sunt prezentate câteva metode numerice de rezolvare a ecuaţiilor cu derivate parţiale, iar în utimul capitol sunt prezentate diferite metode pentru integrarea numerică.

CAPITOLUL 1Erorile de calcul numeric

5

1.1. Surse de eroriSuntem în posesia unui număr suficient de mare de metode numerice

pentru a considera mai în detaliu problema erorilor de calcul numeric. Se observă că o formulă de calcul numeric se aplică de obicei în mod repetat. În consecinţă, prezintă importanţă nu numai eroarea introdusă într-o etapă, ci şi tendinţa de a amplifica sau, dimpotrivă, de a atenua erorile introduse anterior, adică stabilitatea metodei numerice.

Erorile inerente sunt erorile legate de cunoaşterea aproximativă a unor valori provenite din măsurători sau din faptul că avem de-a face cu numere iraţionale. Rezultatul oricăror calcule depinde şi de precizia datelor introduse iniţial. Ca erori inerente pot fi considerate şi erorile de conversie făcute la trecerea în baza 2 a unor numere care se introduc în memoria calculatoarelor numerice. De exemplu, numărul 0.1 reprezentat printr-un număr finit de zecimale în baza 10, devine o fracţie zecimală periodică în baza 2 (0.110 = 0.0(0011)2).

Erorile de metodă sau erorile de trunchiere sunt provenite din aproximaţiile făcute la deducerea formulelor de calcul. De exemplu, restul la interpolarea polinomială, distanţa la rădăcină, din metodele iterative de calcul, etc. Spre deosebire de erorile inerente, erorile de metodă pot fi reduse, în principiu, oricât de mult.

Erorile de rotunjire sunt legate de posibilităţile limitate de reprezentare a numerelor în calculatoarele numerice. În calculator se pot reprezenta numere cu un număr de cifre semnificative în funcţie de lungimea cuvântului (măsurată în biţi) utilizat la stocarea unui număr.

În memoria internă a unui calculator numeric se foloseşte reprezentarea în virgulă mobilă, în formă normalizată. Orice număr real x se scrie

, unde m este un număr real numit mantisă, ( ) este baza sistemului de numeraţie, iar e (întreg) este exponentul.

În forma normalizată, mantisa este cuprinsă în intervalul ().

Excepţia de la această regulă de reprezentare este numărul zero.Deci, un număr real cu mai multe cifre semnificative este „rotunjit” la

numărul de cifre maxim, lucru care se realizează prin rotunjirea mantisei. Alte rotunjiri se fac în decursul operaţiilor.

Notând cu x valoarea exactă a numărului şi cu valoarea calculată (aproximativă), eroarea absolută se defineşte ca

,iar raportul se numeşte eroare relativă, notată cu ,

6

.Fie k numărul de cifre semnificative. Presupunem că b = 10. Atunci, un

număr x se va scrie, ,

unde n conţine cifrele care nu pot incluse în mantisa m. Rotunjirea se face simetric (de obicei), adică se înlocuieşte

dacă , dacă . În acest fel marginea erorii relative este

. Erorile cu marginea dată mai sus se fac la introducerea numerelor reale în memoria calculatorului numeric.

1.2. Propagarea erorilor de calculFie două numere, x şi y, introduse cu erorile , respectiv .

, . Propagarea erorilor la adunare

Presupunem că se efectuează suma numerelor,

astfel încât eroarea relativă la sumare este,

adică o sumă ponderată a erorilor introduse la reprezentarea în calculator a cantităţii sumate. Notăm cu eroarea introdusă suplimentar la reprezentarea sumei . Eroarea relativă totală la sumare, , va fi

.Propagarea erorilor la scădere

Presupunem că se efectuează scăderea numerelor,

astfel încât eroarea relativă la scădere este,

adică o diferenţă ponderată a erorilor introduse la reprezentarea în calculator a diferenţei. Notăm cu eroarea introdusă suplimentar la reprezentarea diferenţei . Eroarea relativă totală la scădere, , va fi

.

Fie , şi . Să se calculeze .

Avem , , deci

, .

7

.

Această eroare relativă a lui se propagă în toate calculele ulterioare.

Propagarea erorilor la înmulţirePresupunem că se efectuează produsul numerelor

,unde s-a neglijat produsul considerat ca având un ordin de mărime suficient de mic. Eroarea la înmulţire este

. Deci la înmulţire erorile relative introduse iniţial se adună. Pot apărea noi erori, deoarece produsul xy poate avea un număr de cifre semnificative mai mare decât cel admis, necesitând o nouă rotunjire. Notând cu această nouă eroare, vom obţine eroarea relativă totală la înmulţirea a două numere

,iar ca margine a erorii

,acoperitoare deoarece erorile se compun după legi probabilistice. Eroarea totală la calculul expresiei a cărei valoare aproximativă este este

,cu mărginirea

. Propagarea erorilor la împărţire

Presupunem că se efectuează împărţirea numerelor

.

Deoarece seria este convergentă, exprimând-o liniar, obţinem

.

Eroare la împărţire este

.

Notând cu noua eroare datorată împărţirii şi cu eroarea totală la împărţirea a două numere obţinem .

8

1.3. Algoritmi şi complexitate de calculProcesele de prelucrare automată a informaţiei se caracterizează prin

natura lor algoritmică şi faptul că sunt efectuate cu ajutorul unor maşini automate de calcul. Calculatoarele electronice sunt capabile să rezolve problemele de calcul descompuse în operaţii elementare, precizându-se succesiunea operaţiilor. Din această perspectivă, prin algoritm se înţelege o mulţime finită de reguli de calcul (care constituie paşii algoritmului), care indică operaţiile elementare necesare rezolvării unei probleme şi ordinea efectuării lor. Orice algoritm porneşte de la datele iniţiale, pe care le prelucrează în vederea obţinerii rezultatului problemei. Algoritmii se caracterizează prin generalitate, finitudine şi unicitate. Oricărei probleme care admite o formulare matematică i se poate asocia un algoritm de rezolvare, care nu se confundă nici cu formularea matematică şi nici cu reprezentarea sub formă unui program într-un limbaj de programare. Printre cele mai utilizate metode de descriere a algoritmilor se numără reprezentarea prin scheme logice şi cea cu ajutorului limbajelor de tip pseudocod. În capitolele următoare vom folosi pentru descrierea algoritmilor un limbaj de tip pseudocod, iar în secţiunile 1.4 şi 1.5 a acestui capitol vom prezenta, pe scurt, câteva din metodele de programare şi câteva din elementele limbajului de programare C în vederea implementării algoritmilor în acest limbaj de programare. Algoritmii prezentaţi sunt exprimaţi într-un pseudocod care posedă următoarele operaţii simple: atribuirea (variabilă expresie), apelul de procedură (funcţie) (execută nume (listă_parametri)), citirea (citeşte listă_variabile), scrierea (scrie listă_expresii) şi în următoarele structuri de control: secvenţa ({instrucţiune_1 ... instrucţiune_n}), selecţia unară (dacă condiţie atunci instrucţiune), selecţia binară (dacă condiţie atunci instrucţiune_1 altfel instrucţiune_2), selecţia multiplă (alege selector dintre valoare1:intrucţiune1 … valoaren:instrucţiunen), iteraţia cu număr prestabilit de ciclări (pentru contor expresie_1:expresie_2 repetă instrucţiune), iteraţia pretestată (cât timp condiţie repetă instrucţiune), iteraţia posttestată (repetă instrucţiune până când condiţie).

Vom înţelege prin problemă algoritmică o funcţie total definită , unde I este mulţimea informaţiilor iniţiale (intrările problemei) iar

F este mulţimea informaţiilor finale. Presupunem că I şi F sunt cel mult numărabile. Dacă este precizat, atunci determinarea lui se numeşte o instanţă a problemei P. Vom folosi pentru o instanţă notaţia p şi prin abuz de notaţie vom scrie „ ”. Un algoritm care rezolvă problema P va porni de la o codificare a unei instanţe oarecare a problemei P şi va oferi o codificare a rezultatului.

Din multitudinea de algoritmi existenţi pentru rezolvarea unei probleme se va alege cel mai „performant”, adică să fie uşor de înţeles,

9

codificat, modificat, depanat şi să utilizeze în mod eficient resursele. Eficienţa unui algoritm este evaluată prin timpul consumat în unitatea centrală şi memoria ocupată, iar pentru algoritmii cu specific numeric mai sunt consideraţi şi factori precum precizia şi stabilitatea numerică. Deoarece analizăm performanţele unui algoritm şi nu performanţele unor calculatoare şi nici software-ul folosit, vom exprima timpul consumat în număr de operaţii (în virgulă mobilă). Vom nota prin numărul de operaţii efectuate de algoritmul A pentru rezolvarea instanţei p şi vom asocia problemei rezolvate un număr numit dimensiunea problemei – reprezentând numărul datelor de intrare.

Resursa timp se exprimă ca o funcţie de dimensiunea problemei, funcţie numită complexitatea în timp a algoritmului . Această funcţie este polinomială sau exponenţială, iar dacă dimensiunea problemei creşte la infinit se obţine complexitatea asimptotică în timp a algoritmului. Comportarea în cazul cel mai nefavorabil a algoritmului A pe o intrare de dimensiune n este

şi ,iar comportarea medie a algoritmului este

şi

(adică media variabilei aleatoare ). Deoarece determinarea exactă a lui este dificilă vom căuta margini superioare şi inferioare pentru . Pentru funcţia ,

. Vom spune că P are complexitatea dacă există un algoritm

A care rezolvă problema P şi are complexitatea . O problemă P are complexitatea dacă orice algoritm A care rezolvă problema P are complexitatea . Un algoritm A cu

, unde p este un polinom în n, se numeşte polinomial. Un algoritm care nu este polinomial se numeşte exponenţial. Astfel algoritmii cu complexitate exponenţială sunt irealizabili, în timp ce algoritmii polinomiali de grad mai mare decât 2 sunt nepractici. Importanţa practică a acestora rezultă şi din următorul exemplu. Presupunem că un pas necesită secunde, adică secunde, atunci pentru n = 40,

un algoritm cu funcţia de complexitate n necesită 0.00004 secunde, un algoritm cu funcţia de complexitate necesită 1.7 minute, un algoritm cu funcţia de complexitate necesită 129 ani (!).

10

1.4. Metode de programareVom prezenta un scurt istoric al unor metode de programare, dar nu

exhaustiv. Am considerat necesar acest lucru, în vederea implementării ulterioare, în limbaje de evoluate de programare, a algoritmilor prezentaţi în lucrare. O clasificare cronologică ar fi următoarea: programarea artizanală, procedurală, modulară, structurată, prin abstractizarea datelor şi orientatǎ pe obiecte.

Programarea artizanală este prima modalitatea de programare, în care iniţiativa şi experienţa programatorului joacă un rol important, fiecare programator îşi are propriile reguli de programare, iar programele au un singur corp de instrucţiuni, sunt lungi şi greu de înţeles.

Programarea procedurală are la bază utilizarea procedurilor (care se mai numesc şi subprograme sau subrutine – unităţi distincte de program) care trebuie parametrizate cu anumite variabile numite parametri formali a căror valori de apel se numesc parametri efectivi. Procedurile realizează o abstractizare prin parametri (ele trebuie să fie generale deci procesarea să facă abstractizare de valorile datelor). La apelare o procedură funcţionează după principiul cutiei negre (se cunosc intrările şi ieşirile, rezultatele din acestea dar nu şi modul de transformare). Sunt proceduri care definesc o valoare de revenire (numite şi funcţii) şi proceduri care nu definesc o valoare de revenire.

Programarea modulară are la bază elaborarea programelor pe module. O colecţie de funcţii înrudite (funcţiile obţinute în urma programării unei subprobleme – componentă a descompunerii arborescente a problemei date în subprobleme mai simple), împreună cu datele pe care le prelucrează în comun formează un modul. Deoarece o parte din datele utilizate în comun de funcţiile modulului nu sunt necesare şi în alte module ele sunt protejate (ascunse în modul).

Programarea structurată. Descompunerea unei probleme în subprobleme mai simple se face succesiv în mai multe etape, până când subproblemele sunt direct programabile sub forma unor proceduri sau module. Această descompunere, numită rafinare pas cu pas, este o descompunere arborescentă. În cadrul procedurilor se folosesc anumite structuri de control a execuţiei. Prelucrarea de bază în cadrul acestor structuri de control este cea de atribuire. Această abordare a programării s-a născut din necesitatea eliminării saltului necondiţionat, determinându-l pe Dijsktra să impună D-structurile de control:

secvenţa (în care prelucrările se execută în ordinea scrierii lor); iteraţia pretestată (cât timp o condiţie este adevărată execută o

prelucrare, prelucrare ce trebuie să modifice valoarea de adevăr a condiţiei);alternativa simplă (dacă o condiţie este adevărată execută o

prelucrare).

11

S-a demonstrat (Bohm şi Jacopini) că orice algoritm se poate descrie doar cu D-structurile, dar pentru o mai bună lizibilitate şi înţelegere a programelor sursă s-au adăugat şi iteraţia posttestată, iteraţia cu un număr prestabilit de ciclări, alternativa binară şi alternativa generalizată.

Programarea prin abstractizarea datelor propune metodologii în care conceptele deduse din analiza problemei să poată fi reflectate cât mai fidel în program. În general un tip abstract de date are două componente: datele membru – care reflectă reprezentarea tipului şi funcţiile membru – care reflectă comportamentul tipului.

Programarea orientată spre obiecte. În cadrul structurilor ierarhice, unde conceptele sunt strâns legate între ele, nu este suficientă doar legătura dintre concepte în care datele membru ale unei clase pot fi obiecte ale unei alte clase. Exprimarea ierarhiilor conduce la atribute suplimentare cum sunt cele de moştenire, atribute care conduc la un model nou de programare numit programare orientată obiect. În vârful ierarhiei se află fenomenul sau forma de existenţă care are trăsături comune pentru toate celelalte componente ale ierarhiei respective. Pe nivelul următor al ierarhiei se află componente care pe lângă trăsăturile comune de pe nivelul superior, mai au trăsături specifice. O ierarhie are, de obicei, mai multe nivele, iar situarea unui element pe un nivel sau altul al ierarhiei este uneori o problemă complexă.

CAPITOLUL 2Rezolvarea numerică a ecuaţiilor şi sistemelor de

ecuaţii algebrice neliniare

2.1. IntroducereForma generală a unei ecuaţii (sistem de ecuaţii) algebrice neliniare

este:

unde f este o funcţie vectorială pe , , cu n componente: , iar x este vectorul necunoscutelor cu n

componente . Pentru , notând , avem o ecuaţie algebrică neliniară sau transcendentă, cu o necunoscută.

2.2. Rezolvarea numerică a ecuaţiilor neliniare Fie o funcţie dată (în anumite cazuri, f va fi o funcţie de la C la C). Presupunând că avem funcţia reală de argument real

,

12

orice valoare pentru care se numeşte rădăcină (soluţie) a ecuaţiei sau zero al funcţiei f (x).

În general, prin rădăcină aproximativă a ecuaţiei se înţelege o valoare apropiată de valoarea exactă . O rădăcină aproximativă poate fi definită în două moduri: numărul cu proprietatea sau, numărul cu proprietatea . Să se găsească una sau mai multe soluţii ale ecuaţiei . Presupunem că f este o funcţie continuă şi cu ajutorul, fie matematic, fie experimental, se va determina un interval [a,b] în care ecuaţia are o rădăcină şi numai una care se va nota cu .

2.2.1. Metoda aproximaţiilor succesiveMetoda constă în construirea unui şir de aproximaţii succesive (şir de

iterare) a lui , rădăcina exactă a ecuaţiei , de forma: = aproximaţia iniţială (dat); ,care să tindă la . Pentru a genera acest şir, se înlocuieşte ecuaţia printr-o ecuaţie echivalentă, în intervalul considerat, , unde g este o funcţie

continuă. De exemplu, sau cu . Apoi, plecând

de la ales în [a,b], se construieşte:

Geometric, se înlocuieşte căutarea intersecţiei graficului funcţiei f cu axa absciselor, prin căutarea intersecţiei curbei de ecuaţie cu dreapta de ecuaţie . Există numeroase moduri de a obţine metoda de aproximaţii succesive, adică de a determina g plecând de la f. În acest sens sunt necesare răspunsuri la întrebările:

i) şirul este convergent?ii) dacă şirul converge, limita sa este ?

Dacă nu-i convergent atunci metoda aleasă trebuie eliminată.

Dacă (unde , deoarece şi cum g este

o funcţie continuă urmează că , adică ) atunci . Deci răspunsul la întrebarea ii) este orice şir.

Deoarece, din punct de vedere numeric se efectuează un număr finit de iteraţii, este necesar să ne preocupe şi următoarele probleme:

13

iii) dacă precizia este dată atunci este necesar să se cunoască cum se opresc iteraţiile astfel încât această condiţie să fie îndeplinită.

iv) deoarece se cere obţinerea rapidă a rezultatului aproximativ, este necesar să fie estimată maniera în care evoluează eroarea în cursul iteraţiilor. Metoda aproximaţiilor succesive conduce la următorul algoritm. Algoritmul 2.2.1.1.

Intrări: g = funcţia din relaţia = precizia

x = aproximaţia iniţialănr_max = numărul maxim de iteraţii admis

Ieşiri: x = soluţia aproximativă îndeplinind sau nr_max{

cât timp şi nr_max {

}

}.

Pentru a detecta timpuriu procesele divergente şi a nu parcurge în mod inutil numărul maxim de iteraţii nr_max, imediat după ce este calculată noua valoare a funcţiei g (care cu semn schimbat înseamnă corecţia rădăcinii) este comparată cu corecţia din pasul anterior (păstrată). Dacă corecţia creşte în valoare absolută procesul este divergent.

Dăm câteva metode clasice a metodei aproximaţiilor succesive.

2.2.2. Metoda LagrangeMetoda Lagrange se mai numeşte şi metoda de părţi proporţionale şi

constă în înlocuirea graficului funcţiei f restrâns la un interval [a,b] prin dreapta determinată de punctele şi . Dreapta are ecuaţia:

.

Dacă este intersecţia dreptei AB cu axa absciselor atunci:

14

.

Presupunând că atunci şi repetând procedeul de mai sus, înlocuind b cu obţinem:

.

În general avem unde funcţia g este definită prin

,

sau

.

2.2.3. Metoda Newton Metoda Newton se mai numeşte şi metoda tangentei şi constă în

înlocuirea graficului funcţiei f restrâns la un interval [a,b] prin tangenta într-un punct al graficului.

Tangenta la curbă în punctul are ecuaţia:.

Dacă este intersecţia acestei tangente cu axa absciselor atunci:

.

Se reîncepe procedeul precedent cu tangenta în punctul , în general se obţine:

unde .

Remarca 2.2.3.1. Metoda lui Newton este bine definită dacă . Este suficient ca să fie un zero simplu al lui f, căci, dacă este

continuă şi „suficient de aproape” de vom avea . Metoda Newton conduce la următorul algoritm. Algoritmul 2.2.3.1.

Intrări: f = funcţia din metodă fd = derivata funcţiei f

= precizianr_max = numărul maxim de iteraţii admis

x = aproximaţia iniţială Ieşiri: x = soluţia aproximativă îndeplinind sau nr_max

{

15

cât timp şi nr_max {

}

}. Dacă pe parcursul procesului iterativ se anulează derivata funcţiei atunci se execută o iteraţie de „salvare” în care corecţia rădăcinii se calculează fără a o împărţi la derivata funcţiei. Această abordare are avantajul că rezolvă şi zerourilor de ordin superior, în care funcţia şi prima ei derivată au zerouri comune.

2.2.4. O teoremă de punct fixÎn exemplele precedente, avem următoarea situaţie:

dat .

Problema care se pune este de a alege o metodă, deci funcţia g astfel încât şirul să fie convergent la . Este necesar ca soluţia ecuaţiei în intervalul [a,b] să fie aceeaşi cu soluţia ecuaţiei în intervalul [a,b].

Are loc următorul rezultat: Teorema 2.2.4.1. Dacă în intervalul [a,b] funcţia g verifică următoarele condiţii:

a)pentru orice are loc ;b)g este o aplicaţie strict contractantă, adică există un număr real L:

astfel încât pentru , avem:

atunci, oricare ar fi şirul definit prin

converge spre unica soluţie a ecuaţiei cu . Demonstraţie.

Dacă şi sunt două soluţii distincte ale ecuaţiei urmează că

cu . Deci . Deoarece obţinem = .

16

Avem .

Deci sau

Deci

Deoarece urmează că pentru .

Deoarece şirul verifică criteriul lui Cauchy urmează că el este convergent spre limita sa şi cum g este funcţie continuă, deoarece este strict contractantă, această limită verifică .

Dacă atunci . Deoarece elementele şirului aparţin intervalului [a,b] (deoarece ) şi converge,

urmează că limita sa aparţine intervalului [a,b].Se poate evalua eroarea făcând p să tindă spre :

.

Se constată că pentru n fixat, eroarea este cu atât mai mică cu cât L se apropie de zero, în timp ce, dacă L este apropiat de 1 eroarea se diminuează.

Dacă g este o funcţie derivabilă atunci o condiţie suficientă astfel încât g să fie strict contractantă este următoarea.

Propoziţia 2.2.4.1. Fie g o funcţie derivabilă în [a,b]. Dacă

verifică atunci g este o aplicaţie strict contractantă în

intervalul [a,b].Utilizând formula creşterilor finite:

cu ,

unde . Condiţia a) din teorema precedentă trebuie să fie îndeplinită.

17

Într-adevăr. Pentru funcţia , unde , ,

, dar ecuaţia nu are punct fix în intervalul considerat,

deoarece .

Intervalul [a,b] în care ipotezele teoremei de punct fix sunt verificate este, în general greu de determinat.

Dacă se poate calcula (estima) atunci are loc următorul rezultat.Propoziţia 2.2.4.2. Fie o soluţie a ecuaţiei cu funcţie

continuă. Dacă atunci există un interval [a,b] conţinând pentru care şirul definit prin şi converge la . Demonstraţie.

Fie . Atunci, din continuitatea lui , există un interval conţinând pentru care

şi .

Funcţia g restrictivă la intervalul [p, q] este strict crescătoare şi verifică

.

Arătăm că dacă atunci . Deoarece , este suficient să arătăm că .

Deoarece urmează că , adică .

Similar, deoarece urmează că .

Deci .Luând , avem verificate ipotezele teoremei de punct fix.

Fie . Atunci există un interval astfel încât

şi .

Se arată şi în acest caz că sunt verificate ipotezele teoremei de punct fix.

Şirul converge la din .Dăm fără demonstraţie următorul rezultat.Propoziţia 2.2.4.3. Fie o soluţie a ecuaţiei . Dacă

este continuă într-o vecinătate a lui şi , , şirul definit de şi nu converge spre .

Dacă atunci şirul poate converge sau diverge. Exemplu. Să se găsească condiţia care asigură existenţa unui interval [a,b], conţinând , în care ipotezele teoremei de punct fix sunt verificate pentru:

18

1) metoda aproximaţiilor succesive;2) metoda lui Lagrange;3) metoda lui Newton.

Rezolvare.1) Fie funcţia g definită prin . Dacă f este derivabilă

atunci . Condiţia devine , adică .

2) Fie , şi funcţia g definită prin

.

Dacă f este derivabilă avem

,

de unde

,

deoarece .Dar

cu ,

de unde.

Deci devine .3) Fie , şi funcţia g definită prin

.

Dacă există atunci .

Deci .Deci există totdeauna un interval în care ipotezele teoremei de punct

fix sunt verificate, în cazul metodei tangentei.Este mai uşoară aplicarea metodei Newton apelând la rezultatul

următor.Teorema 2.2.4.2. Dacă verifică

(1)(2) , (3) ,

19

atunci alegând astfel încât , şirul definit prin şi

converge spre unica soluţie a ecuaţiei în intervalul .Demonstraţie. Condiţiile (1), (2) şi (3) asigură existenţa şi unicitatea unei rădăcini

simple în [a,b] a ecuaţiei .Deoarece

urmează că

cu cuprins între şi .Deci

,adică

. Dacă şi sunt de acelaşi semn pe [a,b] atunci pentru

, , şirul este minorat de plecând de la rangul 1.Dacă şi sunt de semne contrare pe [a,b] atunci

pentru , , şirul este majorat de plecând de la rangul 1.1) Fie pentru . Din ipoteză urmează că

.a) Fie . Atunci şirul este minorat de plecând de la

rangul 1.

Cum urmează că ,

. Deoarece nu se anulează decât pentru , oricare din termenii şirului sunt mai mari decât urmează că este de acelaşi semn cu , deci pozitiv. Deci

20

.

Şirul este descrescător, minorat, deci convergent la .b) Fie . Atunci şirul este majorat de plecând de la

rangul 1.

Cum urmează că ,

. Deoarece nu se anulează decât pentru şi oricare din termenii şirului sunt mai mici decât urmează că este de acelaşi semn cu , deci pozitiv. Deci

.

Şirul este crescător, majorat, deci convergent la . 2) Fie pentru . Demonstraţia este similară cu cea de la punctul 1).

2.2.5. Ordinul metodeiÎn afară de convergenţă, trebuie să ştim cât de „rapidă” este

convergenţa şirului definit prin , adică cum se diminuează eroarea în următoarea iteraţie.

Aceasta ne conduce la următoarea definiţie a ordinului metodei.Definiţia 2.2.5.1. Metoda definită prin este de ordin

p, dacă are limită în mulţimea numerelor reale strict pozitive când n

tinde spre .Explicaţia acestei definiţii.

sau

. Metoda este, deci, de ordinul p, dacă şi numai dacă:

şi .

21

Metoda aproximaţiilor succesive este de ordinul 1 (adică convergenţa este liniară)Într-adevăr. Se cunoaşte, din exemplul 1, că .

Dacă atunci metoda Lagrange este de ordinul 1.Rezultă imediat din exemplul 2.

Metoda lui Newton este de ordinul cel puţin 2.Teorema 2.2.5.1.

Dacă este un zero simplu atunci metoda Newton este de ordinul cel puţin doi. Verificare.

Avem .Deoarece

urmează că .

Dacă este un zero simplu al lui f atunci . Dacă atunci altfel ordinul metodei este mai mare decât 2. Pentru metodele de ordinul doi se mai spune încă: convergenţă este pătratică. Ordinul depinde atât de metoda aleasă cât şi de natura zeroului.

Propoziţia 2.2.5.1. Dacă este un zero al lui f de multiplicitate m,

metoda iterativă definită prin cu este de

ordinul cel puţin doi.

Verificare. În loc să considerăm ecuaţia care admite ca soluţie de multiplicitate m, se poate lua care admite ca soluţie simplă.

Metoda Newton este atunci definită prin:

cu ,

adică

.

2.2.6. Accelerarea convergenţei

22

2.2.6.1. Metoda AitkenRemarca 2.2.6.1.1. În metoda aproximaţiilor succesive şi metoda

Lagrange eroarea este de forma,

unde:

, şi .

Verificare.În metoda aproximaţiilor succesive:

, , ,

, .

,

.

. Din exemplul 1 urmează că . În metoda Lagrange:

.

Din exemplul 2 urmează că , unde cu .

, ,

.

Fie .

.

.

Dacă presupunem atunci

Deci ,

23

adică

.

Deci

Adică, dacă atunci soluţia se obţine numai după două iteraţii succesive. Fie . Este de aşteptat ca şirul definit prin:

să aproximeze mai bine rădăcina decât .Teorema 2.2.6.1.1. Dacă este un şir care converge la având

ordinul de convergenţă unu atunci şirul definit prin

converge mai repede, adică .

Demonstraţie. Deoarece are ordinul de convergenţă unu avem

cu şi . Deci

cu .

Deci .

La fel, . Deci

24

Deci .

Metoda definită de se numeşte metoda lui

Aitken.

2.2.6.2. Metoda SteffensenDacă se aplică procedeul de accelerare a convergenţei al lui Aitken

şirului definit prin avem că este o aproximare mai bună a lui decât . Este, deci, natural să continuăm înlocuind cu , apoi de a calcula şi pentru a aplica accelerarea lui Aitken.

Aceasta revine la a calcula succesiv:

deci

.

Se poate defini în funcţie de plecând de la g:

ceea ce se poate scrie:

.

Aceasta este metoda Steffensen, care se poate defini formal prin funcţia G care dă procesul iterativ:

.

Are loc următorul rezultat: Teorema 2.2.6.2.1. Fie o rădăcină simplă a ecuaţiei .

25

Dacă metoda generată de g este de ordinul unu atunci metoda Steffensen este de ordinul cel puţin doi. Dacă metoda generată de g este de ordinul atunci metoda Steffensen este de ordinul . Demonstraţie.

Pentru , în vecinătatea lui , avem

.

Notând cu avem

.

Dacă ordinul metodei este 1, adică , deoarece

atunci dacă .

Deci metoda Steffensen este de ordinul cel puţin 2.Dacă ordinul metodei este , adică atunci

,

,

,adică metoda Steffensen este de ordinul . Metoda Steffensen furnizează următorul algoritm. Algoritmul 2.2.6.2.1.

Intrări: g = funcţia din metoda neaccelerată = precizia

nr_max = numărul maxim de iteraţii admis x = aproximaţia iniţială Ieşiri: x = soluţia aproximativă îndeplinind sau nr_max

{

cât timp şi nr_max {

26

}

}.

2.2.7. Metoda poziţiei falseMetoda Newton prezintă inconvenientul de a calcula derivata lui f. Se

poate ca această derivată să fie dificil de calculat. Atunci, se înlocuieşte derivata printr-o aproximaţie în funcţie de valori ale lui f.

De exemplu

.

Aceasta este metoda falsei poziţii (corzii). Atunci avem

,

.

Se remarcă că se obţine în funcţie de şi şi nu numai în funcţie de , ca în cazul metodelor precedente.

Dacă se compară cu metoda Lagrange definită prin:

se remarcă că metoda poziţiei false se obţine din metoda Lagrange înlocuind a cu . Se poate, deci, prevedea că ordinul metodei poziţiei false va fi mai bun decât al metodei Lagrange, dar mai puţin bun decât al metodei Newton.

Are loc Teorema 2.2.7.1.

Ordinul metodei poziţiei false este

,

dacă şi . Demonstraţie.

Avem:

27

,

.

Dacă şi atunci

. Atunci:

.Pentru n tinzând spre infinit urmează că , unde

.Ordinul metodei falsei poziţii este p, dacă avem

cu .

Deci: şi ,

de unde

.

Înlocuind în , obţinem:

,

oricare ar fi k şi , deci p trebuie să verifice,

adică

.

2.2.8. Principiul dihotomieiDacă este relativ uşor de a verifica condiţiile în pentru a asigura

convergenţa şirului de iteraţii, este în general greu de a alege - iteraţia

28

iniţială. Procedeul dihotomiei permite să se determine valoarea în vecinătatea lui .

Fie o rădăcină a ecuaţiei cu f continuă.Dacă există un interval închis [a,b] astfel încât să fie

negativ, există cel puţin o rădăcină în intervalul [a,b]. Dacă în plus f este strict monotonă în acest interval atunci această rădăcină este unică.

Se poate atunci considera punctul , mijlocul intervalului [a,b] şi

calcula şi localiza în noul interval având pentru extremităţi

şi punctul a sau b după cum sau este de semn contrar lui .

La fiecare iteraţie se înjumătăţeşte lungimea intervalului. La începutul

iteraţiei n, lungimea intervalului era .

Această metodă furnizează următorul algoritm.Algoritmul 2.2.8.1.

Intrări: a, b = capetele intervalului f = funcţia

= preciziaIeşiri: x = rădăcina localizată {

cât timp şi {

dacă atunci altfel }

}.

2.3. Rezolvarea numerică a sistemelor neliniare

29

Pentru mai multă claritate, ne limităm la sisteme de două ecuaţii cu două necunoscute. Se constată că teoria este identică cu cea a ecuaţiilor înlocuind cu şi valoarea absolută în R cu norma în . Fie şi două aplicaţii în .

Se caută astfel încât

Domeniul D se va alege astfel încât să existe în D o soluţie şi numai una a sistemului.

2.3.1. Metoda aproximaţiilor succesiveFie sistemul echivalent cu sistemul:

Aceasta ne conduce, luând :

Iterând, obţinem

Punând ,

,atunci iteraţia se scrie:

. Drept normă luăm: .

Dăm fără demonstraţie următorul rezultat: Teorema 2.3.1.1. Fie mulţimea închidă . Dacă

1) ;2) există o constantă L: astfel încât:

, atunci şirul definit prin , este convergent şi limita sa aparţine lui D. O condiţie suficientă pentru ca G să fie o aplicaţie contractantă este dată de următoarea teoremă.

Teorema 2.3.1.2. Dacă şi sunt funcţii cu derivate parţiale continui în D, atunci,

,

unde .

30

Demonstraţie.Inegalitatea este echivalentă cu

, adică

. Utilizând formula creşterilor finite avem:

,unde şi . Deci:

unde .

. Deci

Deci:

Metoda aproximaţiilor succesive conduce la următorul algoritm. Algoritmul 2.3.1.1. Intrări: n = număr de necunoscute g = procedura de evaluare

= precizia nr_max = numărul maxim de iteraţii admis

X = aproximaţia iniţialăIeşiri: X = soluţia aproximativă îndeplinind şi nr_max{

execută cât timp şi nr_max {

execută }

31

}.

2.3.2. Metoda NewtonMetoda Newton în cazul unei ecuaţii poate fi interpretată

astfel: plecând de la o valoare aproximând soluţia , se caută o creştere astfel ca .

De fapt, utilizând formula Taylor şi neglijând termenii de ordinul 2, obţinem:

,

de unde şi deci .

Pentru sisteme se procedează la fel. Se dezvoltă fiecare din cele două funcţii în serie Taylor şi se neglijează termenii de ordinul .

Plecând de la , căutăm astfel ca:

Dezvoltăm în serie Taylor şi neglijăm termenii de ordinul 2:

Deci, trebuie să verifice următorul sistem:

Fie matricea . Dacă

determinantul lui (matricea lui Jacobi) este diferit de zero atunci

sistemul are o soluţie şi numai una . Atunci metoda

Newton devine:

Şi în cazul sistemelor, metoda Newton este de ordin cel puţin 2. Pentru

a evita inversarea matricii jacobiene, avem .

Un algoritm furnizat de metoda Newton este:Intrări: n = numărul de necunoscute

f = procedura de evaluare a lui J = procedură de calcul a matricii jacobiene

= precizia nr_max = numărul maxim de iteraţii admis

32

X = aproximaţia iniţialăIeşiri: X = soluţia aproximativă îndeplinind şi nr_max{

repetă

execută execută execută Gauss /* Eliminarea Gauss execută Sist_Tr /*Rezolvarea sistemelor triunghiulare până când sau k > nr_max

}.

2.4. Exerciţii1) Fie ecuaţia . Precizaţi un interval , în

care se găseşte o singură rădăcină reală, apoi determinaţi cu o eroare această rădăcină folosind metodele:

a) metoda lui Newton;b) metoda coardei.

2) Câţi paşi sunt necesari pentru determinarea rădăcinii reale a ecuaţiei , cu o aproximaţie de , prin metodele:

a) metoda lui Newton;b) metoda coardei.

3) Să se rezolve sistemul neliniar

folosind metoda lui Newton cu iteraţia iniţială .

CAPITOLUL 3Rezolvarea sistemelor algebrice liniare

3.1. Introducere

33

Un sistem de ecuaţii algebrice liniare poate fi scris sub forma

,

unde sunt coeficienţi, sunt necunoscutele sistemului, iar sunt termenii liberi. Sunt posibile situaţiile:

(i) dacă atunci trebuie aleşi n – m parametri pentru a obţine o soluţie;

(ii) dacă atunci se caută o soluţie care să minimizeze

;

(iii) pentru m = n, dacă atunci sistemul are o soluţie unică altfel (adică, ) atunci sistemul poate avea o infinitate de soluţii sau poate să nu aibă nici o soluţie. Metodele de rezolvare sunt de două feluri: directe (sau exacte) în care soluţia este obţinută după un număr de operaţii dinainte cunoscut şi metode iterative, care utilizând o aproximaţie iniţială o îmbunătăţeşte de la o etapă la alta.

3.2. Metode directe Ca exemple de metode directe se pot aminti metoda lui Cramer bazată pe calculul determinanţilor, metoda de eliminare a lui Gauss, metoda factorizării directe. Deşi în esenţă foarte simplă, metoda lui Cramer se dovedeşte cu totul ineficientă sub aspectul numărului de operaţii şi nu este în general implementată.

3.2.1. Metoda de eliminare a lui GaussFie un sistem de n ecuaţii liniare cu n necunoscute scris sub forma, unde A este o matrice pătrată, nesingulară, de dimensiune , iar x

şi b sunt vectori coloană de dimensiune n. Metoda constă în a obţine zerouri sub diagonala principală a matricii A succesiv, întâi pe prima coloană, apoi pe coloana a doua ş.a.m.d., pe ultima linie a lui A rămânând doar coeficientul - modificat de operaţiile de eliminare anterioare. Un pas (etapă) k al metodei constă în obţinerea de zerouri sub diagonala principală, în coloana k a matricei A. Etapa este simbolizată prin indice superior atât pentru matricea A cât şi pentru vectorul b. Notăm, iniţial, şi . În urma pasului al fazei eliminării, sistemul are forma echivalentă

34

La pasul k ( ) este eliminată necunoscuta din ultimele n – k ecuaţii şi sistemul este adus la forma

În final se obţine sistemul,

în care matricea este superior triunghiulară. Sistemul se rezolvă prin metoda substituţiei inverse potrivit relaţiilor

Sistemul poate fi transformat într-un sistem superior triunghiular echivalent (adică cu matricea sistemului triunghiulară superior) folosind rezultatele următoarei teoreme: Teorema 3.2.1.1. Dacă , , cu

sunt nesingulare, atunci există , nesingulară şi inferior triunghiulară astfel încât matricea este superior triunghiulară. Demonstraţie. Fie un produs de matrici

elementare de forma , în care este matricea unitate, este

coloana k a lui , iar , în care ultimele elemente sunt neprecizate deocamdată.

35

Înmulţind la stânga matricea sistemului A cu matricea de transformare T obţinem

,sau exprimat printr-o secvenţă de transformări elementare:

în care fiecare operaţie anulează elementele subdiagonale din coloana k ale matricei :

,unde este coloana i a matricei .

.

Componenta i a coloanei transformate este:

Remarcăm faptul că, în timpul transformării, primele k componente nu se modifică, fiind afectate numai elementele subdiagonale.

Din condiţia de anulare a acestora urmează

, pentru ,

determinându-se vectorul , deci şi matricea .Coloanele ( ) se vor modifica astfel:

,

iar pe componente:

, , .

36

Coloanele ( ) au elementele subdiagonale deja anulate ( pentru ) astfel că transformarea le lasă nemodificate:

. Metoda de triangularizare şi substituţia inversă conduc la un

- algoritm. Algoritmul 3.2.1.1.

Intrări: n = dimensiunea sistemului A = matricea sistemului b = vectorul termenilor liberi

Ieşiri: A = matricea sistemului triunghiular b = termenii liberi ai sistemului triunghiular x = vectorul necunoscutelor {

pentru pentru { } pentru {

suma suma + - suma) / } }. În implementarea algoritmului eliminării gaussiene, la fiecare pas k al eliminării se efectuează împărţirea liniei pivot k la elementul diagonal . În practică, este posibil ca acest element să fie egal cu zero, situaţie în care este imposibilă continuarea algoritmului. Erorile de rotunjire în calculul elementelor matricii sunt cu atât mai mici cu cât elementul pivot este mai mare în valoare absolutǎ. Astfel, se efectueazǎ pivotarea parţialǎ pe coloane, ceea ce pentru pasul k al eliminării revine la a găsi elementul maxim

al matricei situat pe coloana k şi liniile , adică şi

permută între ele liniile şi k. Permutarea a două linii sau coloane dintr-o matrice se poate face prin înmulţirea cu o matrice de permutare care diferă de matricea unitate prin: , . Astfel înmulţirea:

are ca efect permutarea liniilor şi k din matricea A, în timp ce: , produce permutarea coloanelor k şi .

37

Căutarea elementului maxim în valoare absolută în submatricea delimitată de ultimele linii şi coloane duce la performanţe mai bune de stabilitate numerică. Astfel se efectuează pivotarea totală.

3.2.2. Factorizarea LU Fiind dată matricea , a unui sistem nesingular ne propunem să descompunem această matrice în produsul matricei inferior triunghiulare L cu matricea superior triunghiulară U:

.

Motivul acestei descompuneri constă în transformarea sistemului dat în două sisteme liniare: şi a căror matrici sunt

triunghiulare şi deci rezolvarea lor este uşoară. Relaţia ce leagă elementele

celor trei matrici poate fi scrisă , .

Sunt posibile mai multe factorizări LU. Matricele L şi U au în total elemente nenule, în timp ce ne furnizează relaţii. Este

necesară particularizarea, de exemplu, a elementelor diagonale a uneia din matricile L sau U. Dacă se consideră elementele diagonale ale matricii U unitare se obţine factorizarea Crout, în timp ce dacă elementele diagonale ale lui L sunt unitare se obţine factorizarea Doolittle.

Considerăm, întâi, . Atunci . Pentru , avem .

Pentru , avem .

La pasul k ( ) avem, ,

.

, ,

.

Fie, acum , . Pentru , avem

38

Pentru , avem .

La pasul k ( ) avem ,…,n, :

.

, :

Factorizarea Crout conduce la următorul - algoritm. Algoritmul 3.2.2.1.

Intrări: n = dimensiunea matricii A = matricea de factorizat

Ieşiri: L, U (folosesc acelaşi spaţiu de memorie ca A) {

pentru { pentru { } pentru {

} } }. Din teorema 3.2.1.1. are loc următorul rezultat: Teorema 3.2.2.1. Dacă toate submatricele principale ale unei matrici A nesingulare de ordinul n sunt nesingulare atunci există o factorizare LU, unde L este triunghiulară inferior şi U este superior triunghiulară.

3.2.3. Factorizarea Cholesky

39

Fie A o matrice simetrică şi pozitiv definită (adică

pentru orice ). Deoarece factorizarea LU strică simetria matricii A, factorizarea Cholesky menţine simetria matricii A, având şi avantajul unui timp de execuţie mai mic pe calculator. Această factorizare constă în descompunerea lui A sub forma , unde L este o matrice inferior triunghiulară, ceea ce revine la

, .

La primul pas avem .

Pentru , avem .

La pasul k ( ) avem ,

.

Pentru , :

.

3.3. Metode iterative Numărul de operaţii pentru metodele directe de rezolvare a sistemelor algebrice liniare este (n fiind numărul de ecuaţii şi necunoscute), în timp de metodele iterative pot conduce la un număr mai mic de operaţii. Principiul general de rezolvare este similar cu cel de la rezolvarea ecuaţiei , în care ecuaţia este scrisă sun forma

,conducând la procesul iterativ

. Fie

şi . Pentru rezolvarea sistemului algebric de ecuaţii liniare

considerăm clasa de metode iterative

,

40

unde şi R sunt parametrii care definesc metoda iterativă.

Pornind de la un element arbitrar se construieşte un şir unde fiecare element reprezintă o aproximaţie a soluţiei sistemului , dacă această soluţie există. Aceasta constituie o metodǎ iterativǎ de rezolvare a sistemului algebric.

Ne interesează condiţiile în care şirul de aproximaţii converge către soluţia sistemului.

Pentru matricea A introducem notaţiile:

,

,

.

3.3.1. Metoda iterativă JacobiDacă , atunci necunoscuta din ecuaţia i este

.

Construim şirul definit prin formulele de recurenţă

, ,

, iar prima aproximaţie este un element din . Relaţia

de mai sus se mai scrie , sau sub

formă matriceală . Raportat la cazul general, în metoda Jacobi avem şi . Şirul definit prin formulele de recurenţă

, , se rescrie prin , ,

41

unde . Vom stabili în ce

condiţii şirul definit mai sus converge la soluţia exactă a sistemului . Fie x soluţia exactă şi eroarea în aproximaţia k : . Sistemul

se mai scrie sau . Deci

sau, trecând la norme şi notând cu obţinem .

Deci, dacă atunci şirul este convergent la soluţia exactă.Ţinând seama de formulele de recurenţă scrise pe componente şi

folosind respectiv norma maximă, norma unu sau norma euclidiană, condiţia de convergenţă devine:

sau

sau .

Din motive de eficienţă verificăm condiţia

considerând şirul convergent începând cu un rang k.Metoda Jacobi conduce la următorul algoritm.Algoritmul 3.3.1.1. Intrări: n = ordinul matricii sistemului

A = matricea sistemului b = vectorul termenilor liberi y = aproximaţia iniţială / precedentă a soluţiei nr_max = numărul maxim admis de iteraţii = precizia

Ieşiri: x = soluţia (aproximaţia cerută) {

cât timp (k < nr_max) şi ( ) pentru { pentru pentru

42

} }.

3.3.2. Metoda iterativă Gauss-SeidelÎn această metodă, construim şirul definit prin

formulele de recurenţă

,

, ,

,

şi un element din . Formulele de recurenţă se pot scrie sub forma

,

sau sub formă matriceală şi .

Astfel şi . Remarcăm faptul că, în această metodă se folosesc noile valori ale componentelor vectorului necunoscutelor imediat ce au fost calculate. Şirul definit prin formulele de recurenţă

, , ,

se rescrie prin , , unde . Condiţia de convergenţă este:

. Deoarece calcularea inversei necesită operaţii, o condiţie simplă de convergenţă a metodei Gauss-Seidel este: pentru fiecare linie a matricii A, elementul diagonal considerat în modul să fie strict mai mare decât suma tuturor modulelor celorlalte elemente ale liniei respective. Notând cu

43

condiţia precedentă se scrie: . Să demonstrăm convergenţa. Sistemul de ecuaţii se mai scrie

, ,

cu valorile exacte ale necunoscutelor şi nenuli.

Obţinem ,

, sau, trecând la module:

, .

Utilizând norma maximă a vectorului eroare ,

obţinem , sau

cu .

Metoda Gauss-Seidel furnizează următorul algoritm.Algoritmul 3.3.2.1. Intrări: n = ordinul matricii sistemului

A = matricea sistemului b = vectorul termenilor liberi y = aproximaţia iniţială / precedentă a soluţiei nr_max = numărul maxim admis de iteraţii = precizia

Ieşiri: x = soluţia (aproximaţia cerută) {

cât timp (k > nr_max) şi ( ) pentru { pentru pentru

44

} }.

3.3.3. Metoda relaxǎriiFie q R*. Metoda relaxǎrii este datǎ de

,

adicǎ , .

Dacǎ obţinem metoda Gauss-Seidel.Existǎ şi alte condiţii în care procedeele iterative sunt convergente.

Dacǎ şi A este o matrice simetricǎ şi strict pozitivǎ atunci şirul de aproximaţii construit prin metoda relaxǎrii converge cǎtre soluţia exactǎ a sistemului Ax = b.

3.4. Exerciţii1) Să se rezolve sistemele:

a)

utilizând metoda Jacobi;b)

utilizând metoda Gauss;c) utilizând metoda Gauss-Seidel.

2) Analizaţi comparat metodele Jacobi, Gauss-Seidel şi relaxării pentru

, , .

CAPITOLUL 4Rezolvarea numerică a problemelor algebrice de

valori şi vectori proprii

4.1. Abordarea elementară a subiectului

45

La probleme algebrice de valori şi vectori proprii conduc numeroase probleme tehnico-ştiinţifice.

Vom arăta în cele ce urmează cum se ajunge la o problemă algebrică de valori şi vectori proprii, extinzând problema rezolvării unui sistem omogen de ecuaţii.

Astfel luând sistemul:

constatăm că acesta se poate scrie sub forma , unde A este matricea formată din elementele cu şi , iar x este vectorul

coloană . Sistemul nu are alte soluţii decât cea banală exceptând cazul

.Să presupunem că matricea A depinde de un parametru şi anume să

luăm , unde I este matricea unitate.În acest fel putem căuta acele valori pentru care şi apoi

putem căuta soluţiile nebanale ale ecuaţiei .Problema pe care o avem de rezolvat se poate scrie sun forma:

. Deci o problemă algebrică de valori şi vectori proprii se exprimă prin relaţia:

unde A este o matrice pătratică de ordin n, x un vector (nenul) care se numeşte vector propriu (la dreapta) ce trebuie calculat iar un număr care se numeşte valoare proprie ce trebuie calculată.

Valorile proprii pentru problema sunt rădăcinile ecuaţiei (numită ecuaţie caracteristică):

, fiind un polinom de grad n în obţinut din dezvoltarea

determinantului . Având o valoare proprie , un vector propriu corespunzător x este soluţia nebanală a sistemului omogen de ecuaţii

.4.2. Aspecte teoretice generale4.2.1. Clase speciale de matrici

Matricea adjunctă (notată ) este matricea transpusă şi conjugată ( dacă ).

46

Dacă A este o matrice pătratică şi , matricea A se numeşte hermetiană (autoadjunctă).

Se demonstrează că:

O matrice A este unitară dacă . Dacă A reală atunci (transpusa) şi dacă atunci A se numeşte ortogonală. Într-o matrice unitară (ortogonală), coloanele formează o bază ortonormală.

O matrice de rotaţie este o matrice care diferă de matricea unitate numai prin liniile p şi q, unde singurele elemente nenule sunt:

,unde c, s sunt două numere complexe astfel ca . (Se poate lua

, cu . Pentru matrici reale, , ).

Se verifică uşor că , adică matricea de rotaţie este unitară (ortogonală). Norma euclidiană a unei matrici este definită prin

.

Prin definiţie, urma matricei A este

.

Dacă H este nesingulară, matricele A şi se numesc asemenea (similare).

Ecuaţiile caracteristice a două matrici asemenea sunt identice, deoarece:

Dacă x este vector propriu al matricei A atunci Hx este vector propriu al matricii , deoarece:

din .

47

O matrice asemenea cu o matrice diagonală se numeşte regulată.Dacă A este regulată atunci ea admite n vectori proprii liniar

independenţi (o bază de vectori proprii).Dacă şi rezultă . Notăm

. Avem

. Obţinem . Teorema 4.2.1.1. Dacă valorile proprii ale unei matrici sunt distincte, atunci matricea este regulată.

Teorema 4.2.1.2. Pentru orice matrice pătratică complexă A există o matrice unitară P, astfel încât

să fie superior triunghiulară. O matrice A este normală dacă

.Lema 4.2.1.1. O matrice normală şi triunghiulară este matrice

diagonală.Demonstraţie. Fie superior triunghiulară (

). Scriem egalitatea ce rezultă din faptul că A este normală:

.

Deci:

48

pentru ,

pentru ,

…,deci T este diagonală.

Teorema 4.2.1.3. O matrice normală este regulată şi admite o bază ortonormală de vectori proprii şi reciproc.

Demonstraţie. Fie A normală. Din teorema 4.2.1.2. există P unitară astfel încât să fie superior triunghiulară.

Avem

( pentru că P este unitară).

( pentru că A este normală).Deci , adică B este normală. Conform lemei 4.2.1.1. rezultă

că B este diagonală. Cum (P unitară) şi B diagonală rezultă că A este regulată. Coloanele lui P sunt vectori proprii şi formează o bază ortonormală.

Reciproc. Dacă A admite o bază ortonormală de vectori proprii avem , unde P este matricea unitară a acestor vectori proprii şi D este

matricea diagonală a valorilor proprii. Avem (P unitară). Avem

Deci , adică A este normală.

4.2.2. Punerea corectă a problemei

49

Ne punem problema efectului perturbaţiei matricei asupra valorilor proprii ale ei (adică la mici variaţii ale datelor de intrare să vedem efectul asupra datelor de ieşire din metodă).

Teorema 4.2.2.1. Fie A o matrice regulată şi fie H o matrice ale cărei coloane sunt vectori proprii liniar independenţi ai matricei A. Atunci, pentru fiecare valoare proprie a matricei perturbate , unde B are elementele de modul subunitare, există o valoare proprie a matricei A astfel încât

. Demonstraţie. Din A regulată rezultă

. Scriem ecuaţia caracteristică a lui şi rezultă:

. Deci

. Dacă atunci teorema este demonstrată. Presupunem că pentru . Atunci . Deci

este singulară (avem teorema: dacă atunci este nesingulară).

Deci . Utilizând norma spectrală:

.

Deci .

Dar .

Deci

.

4.3. Metoda lui JacobiCazul matricelor complexe hermitiene se reduce la cel al matricelor

reale simetrice.Fie , B, C reale. Fie valoare proprie (în acest caz este

reală) căreia îi corespunde un vector propriu de forma cu .Deci , sau

50

adică

.

Deci problema devine reală dar de dimensiune dublă.A hermitiană şi

. Deci matricea de ordin este reală şi simetrică. Ne limităm la o matrice A reală şi simetrică. Deci

A normală. Conform lemei 4.2.1.3., A este regulată şi admite o bază ortonormală de vectori proprii. A regulată, este asemenea cu o matrice diagonală ce are pe diagonală valorile proprii.

Metoda lui Jacobi utilizează matrice de rotaţie reale şi construieşte un şir de matrice asemenea cu A.

Luăm şi în general , , în care

o matrice de rotaţie reală, având la intersecţiile liniilor şi

cu coloanele şi elementele , şi respectiv , , iar în rest elementele matricei unitate. Pentru un k fixat notăm . Alegem p şi q astfel încât

(maximul în modul dintre elementele nediagonale).

Alegem astfel încât să se anuleze pentru elementele simetrice cele mai mari în modul.

Dacă este simetrică atunci

. Prin inducţie

simetrică.

51

Se observă că .

Deci asemenea cu .

Facem calculele în relaţia

şi avem

Deci

celelalte elemente rămânând neschimbate: ,

pentru toţi .Din condiţia de anulare a lui obţinem:

52

.

Alegem cu .

Lucrăm cu numărătorul şi numitorul separat (pentru a evita pierderea preciziei la împărţire).

şi , unde c = numărătorul cu semnul fracţiei, d = numitorul în modul. Teorema 4.3.1. Şirul matricelor construit cu metoda lui Jacobi este convergent şi

unde sunt valorile proprii ale matricii A.Demonstraţie. Să luăm

,

unde este matricea simetrică a elementelor nediagonale ale lui cu zero pe diagonală. Din relaţiile de mai sus rezultă

adică suma pătratelor elementelor nediagonale, cu excepţia celor din poziţiile (p,q) şi (q, p) nu se modifică.

Dar elementele din poziţiile (p,q) şi (q, p) în sunt zero.Deci

pentru că . Rezultă că

cu .

Într-adevăr, pentru . În sunt termeni, deci

53

Deci cu . Rezultă . Deci la limită se

obţine o matrice diagonală. Deoarece toate matricile şirului sunt asemenea cu A elementele diagonale sunt valorile proprii ale lui A. Convergenţa este cel puţin comparabilă cu cea a lui , din punct de vedere al rapidităţii.

Dacă ultima rotaţie care se efectuează până elementele nediagonale devin inferioare preciziei a calculelor ( fiind dat) este atunci avem:

. Cu această precizie vectorii proprii sunt coloanele matricii

.Iniţializăm V = I şi la fiecare pas înmulţim la stânga pe V cu rotaţia

respectivă. Notăm . Înmulţirea de mai sus se transpune în formulele:

Metoda Jacobi conduce la următorul algoritm.Algoritmul 4.3.1.

Intrări: n = ordinul matricei simetrice A A = matricea simetrică = precizia

Ieşiri: x = vectorul cu valorile proprii B = matricea ale cărei coloane conţine vectorii proprii {

cât timp {

pentru pentru dacă atunci

dacă atunci{

54

} dacă atunci

{ } altfel {

}

pentru { dacă şi atunci

{

pentru dacă şi atunci

{

} } } pentru pentru pentru

55

pentru dacă atunci

} pentru pentru { \* vectorul propriu j pentru }}.

4.4. Probleme de valori proprii generalizateÎn unele probleme ale fizicii matematice intervin ecuaţii de valori

proprii generalizate,

unde atât A cât şi B sunt matrici pătrate.O posibilitate de transformare a ecuaţiei de mai sus într-o problemă

standard, de forma ,

constă în înmulţirea ecuaţiei la stânga cu , presupunând că matricea B este nesingulară:

. Dacă A şi B sunt simetrice, nu întotdeauna este simetrică.

Dacă A şi B sunt simetrice şi B este pozitiv definită, se poate transforma problema generalizată într-o problemă standard pornind de la descompunerea Chalescki a matricii B:

,unde L este inferior triunghiulară.

Obţinem:,

sau,

,

adică

pentru matricea simetrică

56

.

Valorile proprii ale matricei coincid cu cele ale lui A, iar vectorii proprii corespunzători sunt

,

adică

.

4.5. Exerciţii1) Să se determine valorile proprii şi vectorii proprii pentru matricea

.

2) Să se determine cea mai mare valoare proprie pentru

.

3) Fie matricile:

a) b) .

Să se determine valorile proprii şi vectorii proprii prin metoda Jacobi. De asemenea, să se scrie programul în C pentru metoda mai sus menţionată.

57

CAPITOLUL 5Aproximarea funcţiilor prin polinoame

5.1. IntroducerePentru o funcţie problema aproximării ei printr-un

polinom se pune fie când este dificil de evaluat f, fie când nu se cunoaşte expresia analitică a lui f ci doar valorile ei în anumite puncte ,

, obţinute în general ca urmarea a unor măsurări şi prezentate într-un tabel. Mulţimea de puncte , , cu proprietatea

, o vom nota prin şi o vom numi diviziune a intervalului . Punctele , vor fi numite nodurile diviziunii.

Alegerea unui polinom pentru aproximarea funcţiei f se justifică prin modul simplu computerizat de evaluare a valorii polinomului într-un punct.

Weierstrass, în 1885, demonstrează că „orice funcţie continuă pe un interval este limita uniformă pe a unui şir de polinoame” şi dă un procedeu pentru calculul acestor polinoame. Bernstein, în 1912, dă o metodă prin care putem calcula, pentru o funcţie continuă pe , un şir

de polinoame uniform convergent către f.Fiind dată o funcţie arbitrară se pune problema stabilirii

unui criteriu de alegere a polinomului P care să aproximeze funcţia dată. Astfel este necesar un instrument matematic pentru măsurarea distanţei dintre două funcţii, care impune situarea într-un spaţiu vectorial normat complet şi în funcţie de stabilirea criteriului de alegere a polinomului P care să aproximeze funcţia f vom avea metode de aproximare a funcţiilor prin polinoame şi anume: aproximarea prin interpolare, aproximarea în medie pătratică.

5.2. Aproximarea prin interpolare Fie şi diviziunea

: . Presupunem cunoscute valorile funcţiei f în nodurile , ale diviziunii şi (eventual) derivatele de anumite ordine ale funcţiei f în aceste noduri. Se pune problema determinării unui polinom P care are aceleaşi valori ca şi f în noduri şi (eventual) derivatele de anumite ordine ale lui f şi P în noduri să coincidă. Un astfel de polinom P îl vom numi polinom de interpolare ataşat funcţiei f şi diviziunii .

5.2.1. Interpolarea polinomială Lagrange

58

Fie şi , astfel încât să fie toate distincte. Presupunem că se cunosc valorile funcţiei f în nodurile ,

, . Se pune problema determinării unui polinom P astfel încât

pentru . Teorema 5.2.1.1. Există un polinom şi numai unul de grad cel mult n astfel ca pentru . Demonstraţie. Dacă P şi Q sunt două polinoame de grad cel mult n astfel ca pentru atunci polinomul este de asemenea de grad cel mult n şi pentru . Deoarece polinomul R se anulează de cel puţin (n + 1) ori şi este de grad cel mult n, urmează că R este polinomul identic nul, adică .

Existenţa polinomului P o demonstrăm construind efectiv polinomul P răspunzând condiţiilor pentru . Să găsim un polinom

în spaţiul vectorial al polinoamelor de grad cel mult n, notat , astfel ca pentru , .

Acest polinom se anulează de n ori pentru , . Deci

.Deoarece urmează că

.

Cele (n + 1) polinoame ( ) sunt liniar independente în spaţiul vectorial al polinoamelor de grad cel mult n. Într-adevăr, fie o combinaţie liniară nulă:

.Pentru , , avem

. Deci pentru . Cele (n + 1) polinoame formează deci o bază a acestui spaţiu vectorial de dimensiune n + 1. Polinomul P căutat se obţine ca o combinaţie liniară de polinoame , adică polinomul dat de formula Lagrange:

.

Într-adevăr,

, pentru .

59

Pentru calculul valorii polinomului de interpolare Lagrange într-un punct a folosind relaţia

a fost definit algoritmul următor: Algoritmul 5.2.1.1.

Intrări: n = gradul polinomului de interpolarea = valoarea argumentului polinomuluix = tabloul absciselor punctelor de interpolarey = tabloul ordonatelor punctelor de interpolare

Ieşiri: = valoarea polinomului de interpolare în a{

pentru {

pentru dacă atunci }

}. Complexitatea metodei este , mai exact operaţii ( adăugări, înmulţiri şi împărţiri). În continuare dăm o altă manieră de scriere a polinomului Lagrange. Fie w polinomul de grad definit prin:

. Atunci

. Deci

cu . Înlocuind prin expresia sa în formula Lagrange, obţinem:

.

În cazul când pentru orice k avem relaţia . Cum

60

obţinem formulele baricentrice

.

În acest caz numărul de operaţii este de ordinul , dar o nouă

evaluare a lui , dacă se conservă valorile lui , nu va cere decât

2n operaţii.

5.2.2. Algoritmul AitkenAtunci când se doreşte modificarea punctelor de interpolare, se

preferă următorul algoritm.Lema 5.2.2.1. Fie Q şi R două polinoame de interpolare de grad cel

mult n, construite respectiv pe şi pe astfel încât: pentru . Atunci polinomul de interpolare de grad n + 1, construit pe

astfel încât pentru poate să se obţină plecând de la Q şi R astfel:

.

Demonstraţie. P este un polinom de grad cel mult n + 1.

,

.

În continuare se va construi polinomul P de grad cel mult n astfel încât: pentru .

61

Considerăm următorul tabel triunghiular:

0 1 2 3 j j + 1 n

0

1

2

3

k

n

Fie următoarea relaţie de recurenţă: pentru

; pentru pentru

.

Teorema 5.2.2.1. Polinomul de interpolare este egal cu . Demonstraţie. Se va arăta că cu este polinomul de interpolare pe punctele , ceea ce va implica faptul că

este polinomul căutat.Facem demonstraţia prin recurenţă pe al doilea indice

pentru . Dacă ipoteza este verificată pentru j fixat atunci pentru : este polinomul de interpolare pe punctele ; este polinomul de interpolare pe punctele . Aplicând lema 5.2.2.1., este polinomul de interpolare pe punctele . 5.2.3. Evaluarea restului la interpolarea Lagrange

62

Fie şi fie . Ce se poate spune despre restul , unde P este polinomul de interpolare, astfel că ,

.Se remarcă faptul că x nu poate să aparţină lui . În aceste caz,

se spune că avem o extrapolare.Aproximarea realizată cu ajutorul polinomului de interpolare

Lagrange este nelocală, în sensul că la evaluarea funcţiei de interpolat contribuie informaţii din toate punctele de interpolare, nu numai din cele din vecinătatea argumentului considerat. Aceasta conduce la oscilaţii complet străine de comportarea reală a funcţiei aproximate. Un exemplu în acest sens este pentru aproximarea funcţiei . Distribuirea în mod echidistant a punctelor de interpolare reduce oscilaţiile până aproape la dispariţie.

Teorema 5.2.3.1. Dacă şi atunci

,

unde . Demonstraţie. Fie , definită prin

unde c este o constantă aleasă astfel încât . Deci , unde . Din condiţiile de interpolare , .

Deci F are n + 2 zerouri distincte şi anume: pe . Aplicând teorema Rolle pe acest interval, va avea cel puţin zerouri distincte, cel puţin zerouri distincte. Din aproape în aproape, va avea cel puţin un zero şi deci

. Sub această formă restul nu poate fi determinat în practică, deoarece punctul , care depinde de x şi de nodurile de interpolare nu este cunoscut.

Dacă, însă, se cunoaşte estimaţia , se poate evalua

.

Interpolarea polinomială Lagrange este utilă pentru construirea formulelor de aproximare a integralelor, derivatelor şi ecuaţiilor diferenţiale. În practică, pentru a aproxima o mulţime de puncte de interpolare, se va limita la polinoame de grad cel mult 10.

Se va prefera pentru polinoame de grad mai mare decât 10, să se utilizeze funcţii spline care au cele mai bune proprietăţi de stabilitate.

63

Exemplul 1. Să se construiască polinomul de interpolare Lagrange a funcţiei f definită pe prin .

Rezolvare. Se constată că se aproximează curba la centrul intervalului, iar la extremităţi se produc oscilaţii numite efecte de margine.

Exemplul 2. Să se construiască polinomul de interpolare Lagrange

a funcţiei f definite pe prin .

Rezolvare. Se remarcă faptul că se produc efecte de margine.Exemplul 3. Se dă funcţia f prin următorul tabel:

0 1 21 10

Se cere polinomul Lagrange .Rezolvare.

5.2.4. Diferenţe divizateFormula restului aminteşte oarecum de restul unei alte aproximări

polinomiale a unei funcţii de n ori continuu diferenţiabilă – formula lui Taylor. Vom arăta că şi polinomul lui Lagrange se poate exprima sub o formă similară cu polinomul lui Taylor, în care derivatele se înlocuiesc cu diferenţe divizate.

Diferenţa divizată de ordinul zero în este .

Prin definiţie diferenţa divizată de ordinul întâi pe nodurile distincte şi este numărul notat cu şi este dat de:

.

Diferenţa divizată de ordinul al doilea pe nodurile distincte , şi este:

.

În general, diferenţa divizată de ordinul k pe noduri distincte se defineşte prin recurenţă:

64

în care diferenţele divizate de la numărător sunt de ordin .Rezultatul care urmează precizează că diferenţa divizată depinde de

nodurile în care este definită şi de valorile funcţiei f în aceste noduri.Teorema 5.2.4.1. Pentru orice şi orice noduri distincte

are loc formula:

.

Demonstraţie. Facem inducţie după k. Pentru , formula din enunţ se verifică conform definiţiei diferenţei divizate de ordinul întâi. Presupunem adevărată formula din enunţ pentru k noduri distincte oarecare şi verificăm pentru noduri. Utilizând formula de recurenţă din definiţia diferenţei divizate de ordinul k şi ipoteza de inducţie avem:

Diferenţele divizate definite prin recurenţă furnizează un algoritm de complexitate . Algoritmul 5.2.4.1.

Intrări: n + 1 = numărul de punctex = tabloul celor n + 1 punctey = tabloul valorilor funcţie f în cele n + 1 puncte

Ieşiri: y = diferenţele divizate de ordin 0 .. n în {

pentru

pentru

}.

65

5.2.5. Formula lui Newton de interpolareCu ajutorul diferenţei divizate putem exprima restul la interpolarea

Lagrange sub o altă formă. Avem:

de unde, conform relaţiei din teorema 5.2.4.1. ,

unde (am renotat) . Polinomul Lagrange

se scrie:.

Pentru fiecare m, , polinomul

are gradul cel mult m şi ,

pentru . Deci

,unde

.

Pentru obţinem ; comparând cu expresia lui , pentru şi , deducem

.Se obţine polinomul lui Newton de interpolare:

Calculul polinomului Newton de interpolare se exprimă prin:

Algoritmul 5.2.5.1.Intrări: n = gradul polinomului de interpolare

a = abscisa în care se calculează polinomulx = tabloul punctelor de interpolarey = tabloul valorilor prin f a punctelor de interpolare

Ieşiri: = valoarea polinomului de interpolare în a{

apelează Dif._div. (n,x,y)

66

pentru {

}

}.5.2.6. Diferenţe finite

Nodurile de interpolare sunt echidistante dacă , .

Numim diferenţă finită de ordinul întâi în x .Se verifică direct că este un operator liniar, iar diferenţa finită de

ordinul k în x se defineşte prin relaţia, .

Formula de calcul pentru este următoarea:

.

Pentru demonstraţie, procedăm prin inducţie. Avem:

Făcând în prima sumă schimbarea de indice şi renotând apoi , obţinem:

Prin inducţie după k, se arată că are loc:

.

5.2.7. Formule de interpolare pe noduri echidistante

67

Să presupunem că punctul x în care se face interpolarea pe nodurile echidistante este apropiat de . Făcând transformarea

în polinomul lui Newton de interpolare şi ţinând cont de legătura dintre diferenţe finite şi divizate, obţinem după simplificări:

Aceasta se numeşte formula lui Newton progresivă pe noduri echidistante.

Dacă x este apropiat de , făcând transformarea , , formula lui Newton devine:

Făcând transformarea şi ţinând cont de legătura dintre diferenţe finite şi divizate, obţinem după simplificări:

numit polinomul Newton regresiv pe noduri echidistante.Formula lui Newton progresivă pe noduri echidistante conduce la

următorul algoritm.Algoritmul 5.2.7.1.

Intrări: n + 1 = numărul de puncte h = pasul

x = tabloul punctelor de interpolare = abscisa în care se calculează polinomul

Ieşiri: = valoarea polinomului de interpolare în {

pentru

{

pentru dacă k este par atunci

altfel

pentru

68

{ pentru

dacă este par atunci altfel

} } }.

Exemplu. Se dă funcţia f prin următorul tabel:0 1 2 3

1 10 29 Se cere polinomul Newton progresiv.

Rezolvare. Se întocmeşte un tabel cu diferenţe finite şi se obţine:

.

5.2.8. Interpolarea polinomială HermiteFie şi noduri de interpolare, toate distincte.

Se caută un polinom P astfel încât:, .

Teorema 5.2.8.1. Există un polinom P şi numai unul singur de grad cel mult astfel încât:

pentru . Demonstraţie. Presupunem că există polinoamele P şi Q de grad cel mult astfel ca şi pentru

. Polinomul , care are de asemenea gradul cel mult , admite fiecare ca rădăcină dublă ( ). Deoarece R

are cel puţin rădăcini şi are gradul cel mult urmează că R este polinomul nul şi deci .

Să arătăm că polinomul

,

unde:

cu

69

îndeplineşte condiţiile din enunţ.Deoarece este un polinom de grad n, şi sunt polinoame de grad . Cum

,

,

,pentru avem:

, ,

adică .

,

.Deci:

, ,

pentru avem , , deci . Pentru evaluarea erorii la interpolarea Hermite avem

Teorema 5.2.8.2. Dacă şi atunci

,

unde . Demonstraţia este asemănătoare celei de la interpolare Lagrange.

Polinomul P având expresia

se numeşte polinomul lui Hermite construit pe baza condiţiilor de interpolare, .

Pentru calculul valorii polinomului Hermite într-un punct a a fost definit un - algoritm.

Algoritmul 5.2.8.1. Intrări: n + 1 = numărul nodurilor de interpolare a = abscisa în care se calculează valoarea polinomului

x = tabloul celor n + 1 puncte de interpolare y = tabloul valorilor funcţiei de interpolat în cele n+1 puncte = tabloul valorilor derivatelor funcţiei de interpolat în cele n+1 puncte

70

Ieşiri: = valoarea polinomului Hermite în a{

s pentru

{ pentru

dacă atunci {

s s + +

}

s

}.

5.2.9. Interpolarea prin funcţii splineStudiul interpolării polinomiale ne conduce la două constatări. Prima

este costul calculelor care devine ridicat când gradul polinomului este mare, iar a doua este că se „stăpâneşte” rău efectul de margine dacă intervalul de interpolare este mare.

O altă aproximare este cea locală, adică făcând interpolarea pe subintervale prin funcţii spline.

Fie punctele de interpolare. Pe fiecare subinterval se caută un polinom de interpolare astfel încât condiţiile de

interpolare să fie îndeplinite: pentru .

În plus se cere ca să fie de clasă . De asemenea, pentru :

, continuitatea primei derivate,,

continuitatea derivatei a doua. Condiţiile de interpolare şi cele de clasă impun a căuta de gradul trei. Exprimăm în funcţie de şi

(necunoscute încă).Prin interpolare liniară avem, punând :

71

. Integrând de două ori obţinem:

,unde şi sunt constante de integrare, care se pot calcula ţinând cont de

şi , adică

şi

,

de unde

şi

.

Înlocuind şi prin expresiile lor obţinem:

Punem condiţia ca derivata să fie o funcţie continuă pe intervalul dat.

Avem

,

echivalent cu:

72

,

sau încă pentru :

.

Am obţinut un sistem de ecuaţii liniare cu necunoscute. Precizând valorile lui şi se poate rezolva sistemul liniar tridiagonal.

Exemplu. Pe intervalul fie funcţiile f şi g definite prin

şi . Să se facă interpolarea prin funcţii spline cu 3, 5

şi respectiv 9 puncte echidistante şi să se dea valorile derivatelor la extremităţi.

Rezolvare. Se constată că pentru 9 puncte de interpolare avem o foarte bună aproximare şi nici un efect de margine.

5.3. Aproximarea în sensul celor mai mici pătrate5.3.1. Problema fundamentală a aproximării liniare

Exemplu. Să aproximăm pe pe intervalul printr-un polinom . Se pune problema să se găsească P astfel ca:

a) să fie minimă;

b) să fie minimă;

c) să fie minimă.

Este necesară o măsură numerică a „distanţei” dintre funcţia de aproximat şi polinomul aproximant într-un spaţiu vectorial normat.

Fie X un spaţiu vectorial normat peste corpul K (adesea ). Notăm prin produsul scalar şi prin norma. Fie , n elemente din X liniar independente. Trebuie să aproximăm y printr-o combinaţie liniară de , adică, să găsim

cu , astfel ca să fie cât mai mic posibil.

73

Definiţia 5.3.1.1. Cea mai bună aproximare a lui y printr-o combinaţie liniară de este un element astfel ca:

. Cea mai bună aproximare constă în a găsi minimul normei erorii.

Drept norme se consideră:1) dacă atunci (modulul);2) dacă atunci (norma euclidiană);3) dacă (= spaţiul vectorial al funcţiilor continue pe

) atunci:

a) (norma convergenţei uniforme);

b) (norma convergenţei pătratice).

5.3.2. Teoreme fundamentaleTeorema 5.3.2.1. Fie n elemente liniar independente din

. Pentru orice , problema celei mai bune aproximări are o soluţie.

Verificare. Trebuie arătat că funcţia norma erorii

are un minim. Pentru aceasta, se arată că această funcţie este continuă şi că este suficient a limita mulţimea de variaţie a lui la o mulţime mărginită şi închisă (exteriorul acestei mulţimi nu conţine minimul) pentru a aplica teorema: orice funcţie continuă pe o mulţime închisă şi mărginită din

îşi atinge minimul său.Corolarul 5.3.2.1. Fie , n întreg fixat. Problema: să se găsească astfel ca:

are o soluţie. Această soluţie este unică.

Aceasta se numeşte aproximarea Tchebycheff de grad n a lui f.Corolar 5.3.2.2. Fie , n întreg fixat. Problema: să se găsească astfel ca:

are o soluţie. Această soluţie este unică.

74

Aceasta se numeşte aproximarea continuă în sensul celor mai mici pătrate.

Corolarul 5.3.2.3. Date fiind , puncte distincte

cu , problema: să se găsească astfel ca:

are o soluţie.

Aceasta se numeşte aproximarea discretă în sensul celor mai mici pătrate.

5.3.3. Aproximarea discretă în sensul celor mai mici pătratea) Cazul general Să se determine unde

,ceea ce revine la a determine numere . Cele funcţii liniar independente sunt date. Se pot alege,

( ) sau funcţiile trigonometrice sau cele exponenţiale.

Se cunosc perechi ( ) cu , unde sunt distincte. Fie , .

Să se determine astfel ca să fie minim.

Considerând norma euclidiană în , avem , unde

. Luăm: . , matricea rectangulară cu linii şi coloane.

Avem .Minimul lui R este caracterizat prin , ,

, adică:

Pentru 0 , împărţind prin 2 şi făcând să tindă spre zero obţinem ),(0 UbUay . pentru 0 , împărţind prin 2 şi făcând să tindă spre

zero obţinem ),(0 UbUay . Cele două inegalităţi dau .

75

Notând cu transpusa lui U obţinem , . Deci valoarea minimă a funcţiei R verifică sistemul liniar

, unde este o matrice pătratică de ordin .Acest sistem se numeşte sistemul ecuaţiilor normale.

b) Cazul liniar În acest caz particular, se caută polinomul de grad unu,

astfel ca , . Fie

şi . Sistemul de ecuaţii normale devine:

,

adică . Notând , ,

deoarece , avem

, care, în

statistică, se numeşte covarianţa vectorilor X şi Y şi se notează

.

Dacă , obţinem .

Rezolvarea sistemului dă , .

5.4. Exerciţii1) Să se determine polinomul de interpolare Lagrange care aproximează funcţia f, cu valorile din tabelul alăturat:

0 1 2

76

0 1 27

2) Să se calculeze, folosind polinomul de interpolare Lagrange, funcţia definită de următorul tabel:

0 1,2 2,5 4,0 5,1 6,0 6,5 7,03,00 6,84 14,25 27 39,21 51,00 58,25 66,00

Se va evalua .3) Pentru , , să se determine diferenţele finite pentru

.

4) Fie funcţia definită prin .

a) Construiţi polinomul P de interpolare Lagrange pe punctele 0, 1, 3, 5.b) Calculaţi şi comparaţi cu .c) Calculaţi polinomul Q al lui Hermite astfel ca:

, , , . Deduceţi valoarea lui . Comparaţi cu . Care este concluzia?5) Să se determine valorile polinoamelor de interpolare: Lagrange, Aitken, Newton şi Hermite ce aproximează următoarele funcţii pe nodurile de interpolare corespunzătoare în punctele precizate şi să se calculeze diferenţele finite şi divizate:a) şi nodurile: 0; 2; 4; 6; 8 în 0,1;b) şi nodurile 1; 2; 3; 4 în 1,1;

c) şi nodurile 2; 3; 4; 5; 6 în 4,5.

De asemenea, să se scrie programele în C pentru metodele mai sus menţionate.

CAPITOLUL 6Rezolvarea ecuaţiilor diferenţiale ordinare de ordinul I

6.1. Introducere Ecuaţiile diferenţiale numerice ordinare sunt o parte a analizei numerice care studiază soluţia numerică a ecuaţiilor diferenţiale ordinare (ODE). Această parte este cunoscută de asemenea sub denumirea de integrare numerică, dar unii cercetători rezervă acest termen pentru calculul integralelor. Multe ecuaţii diferenţiale nu pot fi rezolvate analitic, caz în care trebuie să ne mulţumim cu o aproximaţie a soluţiei. Algoritmii studiaţi aici

77

pot fi folosiţi pentru a calcula astfel de aproximaţie. O metodă alternativă este să folosim tehnicile din calcul pentru a obţine o dezvoltare în serie a soluţiei. Ecuaţiile diferenţiale ordinare apar în multe discipline ştiinţifice, de exemplu în mecanică, chimie, biologie şi economie. În plus, unele metode în numerica ecuaţiilor diferenţiale ordinare transformă ecuaţia cu derivate parţiale într-o ecuaţie diferenţială ordinară, care trebuie rezolvată. Problema Vrem să aproximăm soluţia ecuaţiei diferenţiale:

(1)

unde f este o funcţie care duce [t0,∞) × Rd în Rd, şi condiţia iniţială y0   Rd

este un vector dat. Formularea de mai sus se numeşte problema valorii iniţiale (IVP). Teorema Picard – Lindelöf afirmă că există o soluţie unică, dacă f este Lipschitz continuă. În contrast, problemele valorii mărginite (BVP) menţionează soluţia y în mai mult de un punct. Diferite metode trebuie să fie folosite pentru a rezolva BVP, de exemplu shooting method, multiple shooting sau metode globale metoda diferenţelor finite sau metode sintagmă. Ne limităm la ecuaţii diferenţiale de ordinul I (adică, în ecuaţie apare doar prima derivată a lui y). Oricât, o ecuaţie de ordin superior poate fi uşor convertită într-o ecuaţie de ordinul I introducând noi variabile. De exemplu, ecuaţia de ordinul II

poate fi rescrisă ca ecuaţiile de ordinul I:

şi .

6.2. Metode Prezentăm două metode elementare.

6.2.1. Metoda EulerPlecând de la ecuaţia diferenţială (1), înlocuim derivata prin diferenţa finită

(2)

care conduce la formula următoare: (3)

78

Această formulă este de obicei aplicată în felul următor. Alegem un pas h şi construim şirul t0, t1 = t0 + h, t2 = t0 + 2h, … Notăm prin yn o estimare numerică a soluţiei exacte y(tn). motivaţi de (3), calculăm aceste estimări prin următoarea schemă recursivă:

(4) Aceasta este metoda Euler (sau metoda Euler forward, în contrast cu metoda Euler backward, care va fi descrisă în continuare). Metoda este numită după Leonhard Euler care a descris-o în 1768. Metoda Euler este un exemplu de metodă explicită. Aceasta înseamnă că noua valoare este definită în funcţie de ceva ce deja se cunoaşte,

precum .

6.2.2. Metoda Euler backward Dacă, în locul lui (2), folosim aproximarea:

(5)

obţinem metoda Euler backward:

(6)

Metoda Euler backward este o metodă implicită, adică trebuie să rezolvăm o ecuaţie pentru a găsi . Adesea folosim o iteraţie de punct fix sau metoda Newton – Raphson (sau modificată) pentru a obţine asta. Bineînţeles, ia timp să rezolvăm această ecuaţie; acest cost trebuie luat în considerare când selectăm metoda pe care o folosim. Avantajul metodelor implicite precum (6) este că ele sunt de obicei mai stabile pentru a rezolva o ecuaţie mare, adică poate fi folosit un pas h mai mare.

6.3. Generalizări

79

Metoda Euler este adesea nu destul de precisă. Mai precis, ea are doar ordinul 1, acest lucru făcându-i pe matematicieni să caute metode de ordin mai mare. O posibilitate este folosirea nu doar a valorii anterior calculată

pentru a determina , dar şi determinarea faptului ca soluţia să depindă de mai multe valori precedente. Aceasta conduce la aşa – numita metoda multi – pas. Probabil cea mai simplă este metoda Leapfrog care este de ordinul 2 şi, vorbind grosolan, se bazează pe două valori de timp. Aproape toate metodele practice multi – pas sunt din familia metodelor multi – pas liniare, care au forma:

O altă posibilitate este să folosim mai multe puncte în intervalul [tn,tn+1]. Aceasta conduce la familia metodelor Runge – Kutta. Ambele idei pot fi combinate, metodele rezultante fiind numite metode liniare generale.

6.3.1. Metode Runge – Kutta În analiza numerică, metodele Runge – Kutta sunt o familie importantă de metode iterative implicite şi explicite pentru aproximarea soluţiilor unor ecuaţii diferenţiale ordinare. Aceste tehnici au fost dezvoltate în jurul anului 1900 de matematicienii germani C. Runge şi M. W. Kutta.

6.3.1.1. Metoda clasică Runge – Kutta de ordin 4 Fie problema . Metoda clasică Runge – Kutta de ordin 4 (RK4) pentru această problemă este dată de ecuaţiile următoare:

unde este aproximarea RK4 pentru şi

80

Astfel, valoarea următoare ( ) este determinată de valoarea curentă

( ) plus produsul dintre mărimea intervalului (h) şi o pantă estimată. Panta este o medie ponderată a pantelor:

este panta de la începutul intervalului;

este panta de la mijlocul intervalului, folosind panta

pentru a determina valoarea lui y în punctul folosind metoda lui

Euler; este tot panta de la mijlocul intervalului, dar folosind

panta pentru a determina valoarea lui y;

este panta de la sfârşitul intervalului, cu valoarea lui y

determinată folosind . Atunci,

.

Metoda RK4 este o metodă de ordin 4, adică eroarea pe pas este de ordinul lui , în timp ce eroarea totală acumulată are ordinul .

6.3.1.2. Metode Runge – Kutta explicite Familia metodelor Runge – Kutta explicite este o generalizare a metodei RK4. Ea este dată de:

,

unde

Metoda Runge – Kutta este consistentă dacă

pentru .

81

6.3.1.3. Metode Runge – Kutta ajustativeMetodele ajustative produc o estimare a erorii de trunchiere locale a

unui pas simplu Runge – Kutta. Acest lucru se produce având două metode, una de ordin p şi una de ordin .

Pasul cu ordinul cel mai mic este dat de

,

unde sunt aceeaşi pentru metoda de ordinul cel mai înalt. Atunci, eroarea este:

care este .

6.3.1.4. Metode Runge – Kutta implicite Metodelor implicite sunt mai generale decât cele explicite. Soluţia aproximativă a problemei valorii iniţiale reflectă numărul mai mare de coeficienţi:

,

unde

. Cel mai simplu exemplu de metodă Runge – Kutta implicită este metoda Euler backward:

.

6.3.2. Caracteristici O bună implementare a uneia din aceste metode de rezolvare a ecuaţiilor diferenţiale ordinare este necesară mai mult decât formula time – stepping. Adesea este insuficient să folosim acelaşi pas tot timpul, deci au fost dezvoltate metodele care variază pasul. De obicei, pasul este ales astfel încât eroarea (locală) pe pas să fie mai mică decât un nivel de toleranţă. Aceasta înseamnă că metodele trebuie de asemenea să calculeze un indicator al erorii, o estimare a erorii locale. O extindere a acestei idei este să alegem dinamic între diferite metode de diferite ordine (aceasta se numeşte metoda variaţiei ordinului). Metodele

82

de extrapolare sunt des folosite pentru a construi diferite metode de diferite ordine. Alte caracteristici sunt:

ieşirea densă: aproximări numerice mici pentru întreg intervalul de integrare, şi nu doar în punctele t0, t1, t2, ...

locaţia evenimentului: găsirea timpul în care, o funcţie particulară se anulează.

ajutor pentru calcul paralel. când sunt folosite pentru integrare în raport cu timpul, reversibilitatea

timpului.

6.4. Metode alternativeMulte metode nu cad în contextul discutat aici. Unele clase de metode

alternative sunt: metodele multi – derivate, care folosesc nu numai funcţia f ci şi

derivatele ei. Această clasă include metodele Hermite – Obreschkoff şi metodele Fehlberg, precum şi metode precum metoda Parker – Sochacki, care calculează recursiv coeficienţii seriei Taylor a soluţiei y.

metode pentru ODE de ordinul 2. Spunem că toate ODE de ordin superior pot fi transformate în ODE de primul ordin de tipul (1). În timp ce aceasta este cu siguranţă adevărat, s-ar putea să nu fie cea mai bună cale de procedat. În particular, metodele Nyström lucrează direct cu ecuaţii de ordinul doi.

metode geometrice intrinseci sunt special destinate pentru clase speciale de ODE (de exemplu, integratorii simplectici pentru soluţia ecuaţiilor hamiltoniene). Ele ţin cont că soluţia numerică respectă geometria acestor clase.

Analiza Analiza numerică nu este doar schiţa pentru metodele numerice, dar şi pentru analiza lor. Cele trei concepte centrale în această analiză sunt: convergenţa, ordinul şi stabilitatea.

Convergenţa O metodă numerică este convergentă dacă soluţia numerică se apropie de soluţia exactă când pasul h tinde către 0. mai precis, cerem ca pentru fiecare ecuaţie diferenţială ordinară (1) cu o funcţie Lipschitz f şi pentru orice t* > 0,

.

Toate metodele menţionate mai sus sunt convergente. De fapt, convergenţa este o condiţie sine qua non pentru orice schemă numerică.

83

Consistenţa şi ordinul Presupunem că metoda numerică este:

. Eroarea locală a metodei este eroarea făcută de un pas al metodei. Adică, ea este diferenţa dintre rezultatul dat de metodă, presupunând că nu a fost făcută nici o eroare în paşii anteriori, şi soluţia exactă:

.

Metoda este consistentă dacă .

Metoda este de ordinul p dacă pentru .

Deci, o metodă este consistentă dacă are un ordin mai mare decât 0. Metoda Euler (forward) (4) şi metoda Euler backward (6) introduse mai sus au amândouă ordinul 1, deci sunt consistente. Majoritatea metodelor folosite în practică ating un ordin mare. Consistenţa este o condiţie necesară pentru convergenţă, dar nu şi suficientă; pentru ca o metodă să fie convergentă trebuie să fie atât consistentă cât şi stabilă în zero. Un concept asociat este eroarea globală, adică eroarea făcută în toţi paşii pe care trebuie să-i atingem la timpul fixat t. explicit, eroarea globală la

timpul t este unde . Eroarea globală a unei metode cu

un pas de ordinul p este ; în particular, o astfel de metodă este

convergentă. Această afirmaţie nu este neapărat adevărat pentru metodele multi – pas. Stabilitatea şi stiffness Pentru unele ecuaţii diferenţiale, aplicarea unor metode standard – precum metoda Euler, metodele explicite Runge – Kutta sau metodele multi – pas (de exemplu, metodele Adams – Bashforth) – expun instabilitatea în soluţie, totuşi alte metode pot conduce la soluţii stabile. Acest “comportament dificil” în ecuaţie (care nu trebuie neapărat să fie dificilă) este descris ca rigiditate, şi adesea este cauzat la prezenţa unor scale de timp diferite în problema de bază. Problemele rigide sunt omniprezente în cinematica chimică, teoria controlului, mecanica solidului, predicţia vremii, biologie şi electronică.

84

CAPITOLUL 7Integrarea numerică

7.1. Introducere În analiza numerică, integrarea numerică constituie o familie vastă de algoritmi pentru calculul valorii numerice a unei integrale definite, şi prin dezvoltare, termenul este de asemenea numit uneori folosit să descrie soluţia numerică a unei ecuaţii diferenţiale. Termenul cvadratură numerică (pe scurt, adesea, cvadratură) este mai mult sau mai puţin sinonim pentru integrarea numerică, special ca aplicat integralelor unu – dimensionale. Integrarea 2 – dimensională şi mai mare este uneori descrisă ca cubatură, deşi sensul cvadraturii este înţeles la fel de bine şi pentru integrări de dimensiune mai mare. Problema fundamentală considerată de integrarea numerică este să calculăm o soluţie aproximativă a unei integrale definite:

.

Dacă este o funcţie netedă peste o un număr mic de dimensiuni şi limitele de integrare sunt mărginite, atunci există multe metode de aproximare a integralei cu o precizie arbitrară.

Legătura cu ecuaţiile diferenţiale

Problema evaluării integralei poate fi redusă la problema

valorii iniţiale pentru o ecuaţie diferenţială ordinară. Dacă integrala de mai sus se notează cu , atunci funcţia I satisface

. Metode dezvoltate pentru ecuaţii diferenţiale ordinare, precum metodele Runge – Kutta, pot fi aplicate pentru problema formulată şi astfel pot fi folosite pentru a evalua integrala. Această ecuaţie diferenţială are o formă specială: membrul drept conţine doar variabila dependentă (x) şi nu variabila independentă (I). Acest fapt sugerează că pot fi dezvoltate metode specifice pentru evaluarea unei integrale, şi că aceste metode pot lucra mai bine decât metodele generale pentru problema valorii iniţiale pentru ecuaţii diferenţiale. Motive pentru integrarea numerică Există câteva motive pentru a face integrarea numerică. Integrantul

poate fi cunoscut doar în anumite puncte, aşa cum este obţinut prin prelevare. Unele sisteme încastrate şi alte aplicaţii pe calculatoare pot avea nevoie de integrare numerică pentru acest motiv. Poate fi ştiută o formulă

85

pentru integrant, dar este dificil sau imposibil de găsit o antiderivată care este

o funcţie elementară. Un exemplu de asemenea integrant este ,

antiderivata căruia nu poate fi scrisă într-o formă elementară. Poate fi găsită o antiderivată simbolică, dar este mai uşor de calculat o aproximare numerică decât de calculat antiderivata. Acesta poate fi cazul dacă antiderivata este dată ca o serie infinită sau ca product, sau dacă evaluarea ei necesită o funcţie specială ca nu este disponibilă. Metode pentru integrale de dimensiune unu Metodele de integrare numerică pot fi descrise ca, combinând evaluările integrantului pentru a obţine o aproximare a integralei. O importantă parte a analizei oricărei metodei de integrare numerică este studiul comportării erorii aproximării ca o funcţie de numărul de evaluări ale integrantului. O metodă care produce o eroare mică pentru un număr mic de evaluări este de obicei considerată ca superioară. Reducând numărul de evaluări ale integrantului reducem numărul de operaţii aritmetice implicate, şi de aceea, reduce eroarea totală rotunjită. De asemenea, fiecare evaluare ia timp, şi integrantul poate fi complicat. Poate fi făcută integrarea numerică brută dacă integrantul se comportă rezonabil (adică, este continuu) prin evaluarea integrantului cu creşteri foarte mici. Regulile cvadraturii bazate pe interpolarea funcţiilor O clasă largă de reguli de integrare poate fi obţinută prin construcţia funcţiilor de interpolare car sunt uşor de integrat. De obicei aceste funcţii de interpolare sunt polinoame.

Exemplificare pentru regula dreptunghiului Cea mai simplă metodă de acest tip este să punem ca funcţie de interpolare o funcţie constantă (un polinom de grad 0) care trece prin punctul

. Această metodă se numeşte regula mijlocului sau regula

dreptunghiului.

.

86

.Exemplificare pentru regula trapezului

Funcţia de interpolare poate fi o funcţie afină (un polinom de grad 1) care trece prin punctele şi . Această metodă se numeşte regula trapezului.

Exemplificare pentru regula Simpson Pentru oricare din aceste trei reguli, putem face o aproximare mai precisă împărţind intervalul într-un număr de n subintervale, calculând o aproximare pe fiecare subinterval, apoi adunând toate rezultatele. Acest mod se numeşte regulă compusă, regulă extinsă sau regulă iterativă. De exemplu, regula trapezului compusă poate fi exprimată ca:

(*)

unde subintervalele au forma , cu şi .

Interpolare cu polinoame evaluate în puncte echidistante în produce formulele Newton – Cotes, care are ca exemple regula dreptunghiului şi regula trapezului. Regula Simpson, care se bazează pe un polinom de grad 2, este de asemenea o formulă Newton – Cotes. Regula corespunzătoare cu fiecare interval subdivizat include toate punctele curente, deci valorile acelor integranţi pot fi refolosite. Dacă permitem intervalelor dintre punctele de interpolare să varieze, găsim un alt grup de formule de cvadratură, precum formulele de cvadratură gaussiană. O regulă de cvadratură gaussiană este mai precisă decât o regulă Newton – Cotes care cere acelaşi număr de evaluări ale funcţiilor, dacă integrantul este neted (adică, dacă el are multe derivate). Alte metode de cvadratură cu

87

variaţia intervalelor include metoda cvadratura Clenshaw – Curtis şi metoda cvadraturii Fejér. Algoritmi ajustaţi Dacă nu are multe derivate în toate punctele, sau dacă derivatele devin mari, atunci cvadratura gaussiană este adesea insuficientă. În acest caz, un algoritm similar cu următorul va executa mai bine: // Acest algoritm calculează integrala definită a unei funcţii de la 0 la 1, ajustativ, // alegând cei mai mici paşi în vecinătatea punctelor problematice. // Fie initial_h mărimea pasului iniţial. x:=0 h:=initial_h accumulator:=0 WHILE x<1.0 DO IF x+h>1.0 THEN h=1.0-x END IF IF error in quadrature of f(x) over [x,x+h] is too large THEN Make h smaller ELSE accumulator:=accumulator + quadrature of f over [x,x+h] x:=x+h IF error in quadrature of f(x) over [x,x+h] is very small THEN Make h larger END IF END IF END WHILE RETURN accumulator Unele detalii ale algoritmului cer o gândire atentă. Pentru multe cazuri, nu este evidentă estimarea erorii dintr-o cvadratură pe un interval pentru o funcţie . O soluţie populară este folosirea a două reguli diferite de cvadratură şi folosirea diferenţei lor ca o estimare a erorii din cvadratură. Cealaltă problemă este să decidem ce semnifică “prea mare” sau “foarte mic”. Un criteriu local pentru “prea mare” este că eroarea cvadraturii trebuie să nu fie mai mare decât unde t, un număr real, este toleranţa pe care dorim să o fixăm pentru eroarea globală. Apoi, din nou, dacă h este deja mic, s-ar putea să nu merite să facem mai mic chiar dacă eroarea cvadraturii este aparent mare. Un criteriu global este ca suma erorilor pe toate intervalele trebuie să fie mai mică decât t. Acest tip de analiză a erorii se numeşte, de obicei, “a posteriori” din moment ce noi calculăm eroarea după ce am calculat aproximarea.

88

Metode de extrapolare Precizia unei reguli de cvadratură de tip Newton – Cotes este, în general, o funcţie de numărul punctelor de evaluare. Rezultatul este mai precis când numărul punctelor de evaluare creşte, sau, echivalent, când mărimea pasului dintre puncte scade. Este natural să ne întrebăm ce va fi rezultatul dacă permitem pasului să se apropie de zero. La această întrebare se poate răspunde extrapolând rezultatul de la două sau mai multe mărimi diferite de zero (de exemplu, extrapolarea Richardson). Funcţia de extrapolare poate fi o funcţie polinomială sau una raţională. Estimarea (a priori) conservativă a erorii Fie f cu prima derivată mărginită pe . Teorema de medie pentru f, unde x < b, ne dă:

pentru un din depinzând de x. Dacă integrăm în x de la a la b în ambii membri şi luăm valoarea absolută, obţinem:

Putem aproxima, în plus, integrala din membrul drept introducând valoarea absolută în integrant şi înlocuind termenul în printr-o limită superioară:

(**)

Deci, dacă aproximăm integrala prin regula cvadraturii

, atunci eroarea nu este mai mare de membrul drept al (**). Putem transforma aceasta într-o analiză a erorii pentru suma Riemann (*), dând marginea superioară a:

pentru termenul erorii a acelei aproximări particulare. Folosind mai multe derivate şi cvadratura putem face o analiză a erorii similară folosind seriile Taylor (folosind o sumă parţială cu rest) pentru f. Această analiză a erorii dă o limită superioară a erorii, dacă derivatele lui f sunt disponibile. Această metodă de integrare poate fi combinată cu aritmetica intervalului pentru a obţine demonstraţii cu calculatorul şi pentru a verifica calculele.

89

Integrale multidimensionale Regulile cvadraturii discutate până acum sunt destinate pentru a calcula integrale de dimensiune 1. Pentru a calcula integrale în dimensiuni multiple, o abordare este să exprimăm integrala multiplă ca integrale de dimensiune unu repetate aplicând teorema lui Fubini. Această abordare cere ca evaluările funcţiei să crească exponenţial cu numărul cu care creşte dimensiunea. Aceste două metode sunt cunoscute ca învingând această aşa – numită curse a dimensiunii. Monte Carlo Metodele Monte Carlo şi metodele cvasi – Monte Carlo sunt uşor de aplicat integralelor multidimensionale, şi pot produce o mai mare precizie pentru acelaşi număr de evaluări ale funcţiei decât integrări repetare folosind metode unu – dimensionale. O clasă mare de metode Monte Carlo utile este formată din aşa – numiţii algoritmi “lanţul Monte Carlo al lui Marcov”, care include algoritmul Metropolis – Hastings şi prelevarea Gibbs.

Integrarea numerică prin metoda Monte Carlo: nodurile sunt aleatoriu echidistante.

Noile noduri sunt albastru închis, vechile noduri sunt bleu. Valoare integralei tinde către 3.32.

7.2. Integrarea în puncte echidistante7.2.1. Formulele Newton-Cotes În analiza numerică, formulele Newton – Cotes, numite şi regulile Newton – Cotes, reprezintă un grup de formule pentru integrarea numerică (numită şi cvadratură) bazate pe evaluarea integrantului în puncte echidistante. Ele sunt numite după Isaac Newton şi Roger Cotes. Formulele Newton – Cotes pot fi utile dacă este dată valoarea integrantului în punctele

90

echidistante. Dacă este posibil să schimbăm punctele în care este evaluat integrantul, atunci alte metode precum cvadratura gaussiană şi cvadratura Clenshaw – Curtis sunt, probabil, mai potrivite. Este presupus că valoarea unei funcţii f este cunoscută în punctele echidistante xi, pentru . Există două tipuri de formule Newton – Cotes: tipul “închis” care foloseşte valoarea funcţiei în toate punctele, şi tipul “deschis” care nu foloseşte valorile funcţiei în capetele intervalului. Formula Newton – Cotes închisă de grad n este exprimată astfel:

unde xi = h i + x0, cu h (numit pas) egal cu . se numesc ponderi.

Aşa cum se poate vedea în următoarea derivaţie, ponderile derivă din polinoamele Lagrange fundamentale. Aceasta înseamnă că ele depind doar de xi nu şi de funcţia f. Fie polinomul de interpolare în forma Lagrange

pentru punctele . Atunci

Formula Newton – Cotes deschisă de grad n este exprimată ca:

Ponderile sunt găsite într-o manieră similară ca în formula închisă. Instabilitatea pentru grad mare

O formulă Newton – Cotes de orice grad n poate fi construită. Oricum, pentru n mare, o regulă Newton – Cotes poate suferi uneori din cauza fenomenului Runge unde eroarea creşte exponenţial pentru n mare. Metode precum cvadratura gaussiană şi a cvadraturii Clenshaw – Curtis cu puncte ne-echidistanete (grupate în extremităţile intervalului de integrare) sunt stabile şi mult mai precise, şi sunt preferate formulei Newton – Cotes. Dacă aceste metode nu pot fi folosite, deoarece integrantul este dat doar în punctele echidistante, atunci fenomenul Runge poate fi evitat folosind o regulă compusă, aşa cum este explicat dedesubt.

91

Formulele Newton – Cotes închisegrad

Denumire Formulă Eroarea

1Regula

trapezului

2Regula lui Simpson

3 regula lui

Simpson

4Regula lui Boole

sauRegula lui Bode

Exponentul pasului h în eroare arată viteza la care descreşte eroarea aproximării. Derivata lui f în termenul eroare arată care polinoame pot fi integrate exact). De notat că derivata lui f în termenul eroare creşte cu 2 pentru fiecare altă regulă. Numărul ξ se află între a şi b.

Formulele Newton – Cotes deschise

grad Denumire Formulă Eroarea

1 Regula dreptunghiului

2 Fără nume

3 Fără nume

4 Fără nume

92

Reguli compuse Pentru ca regulile Newton – Cotes să fie precise, pasul h trebuie să fie mic, ceea ce înseamnă că intervalul de integrare [a,b] trebuie să fie mic. Din acest motiv se execută integrarea numerică împărţind [a,b] în subintervale mai mici, aplicând o regulă Newton-Cotes pe fiecare subinterval şi adunând rezultatele. Aceasta se numeşte regulă compusă.

7.2.2. Metoda dreptunghiurilor În calculul integral, metoda dreptunghiurilor foloseşte o aproximare a unei integrale definite, făcută prin găsirea ariei unei serii de dreptunghiuri. În calculul numeric, aceasta a fost înlocuită de metode de integrare numerică mai sofisticate. Oricare colţ stâng sau drept, sau mijlocul de sus a dreptunghiului se găsesc pe graficul funcţiei, cu bazele de-a lungul axei x. Aproximarea este luată prin sumarea ariilor celor n dreptunghiuri care umplu spaţiul între două valori x dorite.

,

unde , .

Necesitatea lui a + i'Δx apare când a este nenul, şi cum poziţia primului dreptunghi nu este în f(i'Δx) ci, mai degrabă, f(a + i'Δx). Cu cât creşte n, cu atât aproximarea este precisă. De fapt, limita aproximării pentru

este exact egală cu integrala definită:

. Acest lucru este adevărat indiferent de

ce i este folosit. Deşi aproximarea în mijloc tinde să fie mai precisă pentru n finit, limita celor trei aproximări, cum , este aceeaşi, astfel oricare dintre ele poate fi folosită pentru a calcula o integrală definită.

93

Aproximarea dreaptă Riemann Aproximarea în mijloc Aproximarea stângă Riemann Eroarea Erorii de aproximare în mijloc se descompune ca cubul lăţimii dreptunghiului

pentru oarecare.

7.2.3. Regula trapezului În matematică, regula trapezului este o modalitate de a calcula

aproximativ integrala definită . Regula trapezului lucrează

aproximând regiunea de sub graficul unei funcţii printr-un trapez şi calculând aria lui. Rezultă că

.

Pentru a calcula această integrală precis, împărţim, pentru început, intervalul de integrare în n subintervale mai mici, şi apoi aplicăm regula trapezului pe fiecare dintre ele. Obţinem regula trapezului compusă:

.

Aceasta poate fi scrisă alternativ ca

unde

,

pentru (putem folosi de asemenea o reţea neuniformă). Regula trapezului este una din familiile de formule pentru integrarea numerică numite formulele Newton – Cotes. Regula Simpson este alt membru al aceleiaşi familii, adesea mai precis. Regula Simpson şi alte metode ca ea pot fi aşteptate să îmbunătăţească regula trapezului pentru funcţii care sunt de două ori continue diferenţiabile; oricum, pentru funcţii grosolane, regula trapezului este probabil preferabilă. În plus, regula trapezului tinde să devină extrem de precisă când sunt integrate funcţii

94

periodice pe perioada lor, un lucru cel mai bine înţeles în legătură cu formula de sumare Euler – MacLaurin. Un avantaj al regulii trapezului este acela că semnul erorii aproximării este cunoscut. O integrală aproximată cu această regulă pe o funcţie sus-concavă va fi supraestimată deoarece trapezele includ toată aria de sub curbă şi se extind peste ea. Folosind această metodă la o funcţie jos-concavă produce o subestimare deoarece aria este nenumărată pentru partea de sub curbă, dar nici una nu este numărată deasupra. Dacă intervalul integralei fiind aproximat include un punct de inflexiune, atunci eroarea este mai greu de identificat.

7.2.4. Metoda Romberg În analiza numerică, metoda Romberg (1955) generează o matrice triunghiulară formată din estimări numerice ale unei integrale definite

folosind extrapolarea Richardson în repetate rânduri pe regula trapezului. Metoda Romberg evaluează integrantul în puncte echidistante. Integrantul trebuie să aibă derivate continue. Dacă este posibil să evaluăm integrantul în puncte neechidistante, atunci alte metode precum cvadratura gaussiană şi cvadratura Clenshaw – Curtis sunt mai precise.

Metoda Romberg poate fi definită inductiv astfel:

sau

unde . Eroarea este de ordinul

Prima extrapolare, , este echivalentă cu regula Simpson cu puncte. Când evaluările funcţiei sunt costisitoare, ar fi de preferat să înlocuim interpolarea polinomială a lui Richardson cu interpolarea raţională propusă de Bulirsch şi Stoer (1967). Exemplu

95

Ca exemplu, funcţia Gauss este integrată de la 0 la 1, adică funcţia eroare . Matricea triunghiulară este calculată pe linie şi calculul se termină dacă ultimele două intrări din ultima linie diferă prin mai puţin de 1.e – 8. 0.77174333 0.82526296 0.84310283 0.83836778 0.84273605 0.84271160 0.84161922 0.84270304 0.84270083 0.84270066 0.84243051 0.84270093 0.84270079 0.84270079 0.84270079 Rezultatul din colţul drept inferior al matricei triunghiulare este precis. Este de remarcat că acest rezultat derivă din aproximările mai puţin precise obţinute prin regula trapezului în prima coloană a matricei triunghiulare.7.2.5. Regula Simpson În analiza numerică, regula Simpson este o metodă pentru integrarea numerică, aproximarea numerică a integralelor definite. Anume, ea este

aproximarea .

Interpolarea cvadraticăÎnlocuim integrantul prin polinomul integral P(x) care ia

aceleaşi valori ca în punctele terminale a şi b şi în mijlocul .

Putem folosi interpolarea polinomială Lagrange pentru a găsi o expresie pentru acest polinom,

.

Un calcul simplu arată că

.

Media mijlocului şi regulile trapezului O altă derivaţie construieşte regula Simpson din două aproximări mai simple:

regula mijlocului:

regula trapezului: .

Erorile în aceste aproximări sunt, respectiv,

şi .

96

Rezultă că termenul principal al erorii se anulează dacă luăm media

ponderată . Această medie ponderată este exact regula Simpson.

Folosind altă aproximare (de exemplu, regula trapezului cu un număr dublu de puncte), este posibil să luăm o medie ponderată potrivită şi să eliminăm alt termen al erorii. Aceasta este metoda Romberg.

Coeficienţii nedeterminaţiA treia derivaţie începe de la ansatz

.

Coeficienţii şi pot fi fixaţi cerând ca această aproximare să fie exactă pentru toate polinoamele cvadratice. Aceasta produce regula Simpson. Eroarea Eroarea în aproximarea unei integrale prin regula Simpson este

, unde este un număr oarecare între a şi b. Eroarea este

(asimptotic) proporţională cu . Oricum, derivaţiile de mai sus

sugerează o eroare proporţională cu . Regula Simpson câştigă un

ordin în plus deoarece punctele în care integrantul este evaluat, sunt distribuite simetric în intervalul . Regula Simpson compusă Dacă intervalul de integrare este mic (adică funcţia de integrat este relativ netedă pe intervalul ), atunci regula Simpson va furniza o aproximare adecvată a integralei exacte. Pentru o astfel de funcţie, un interpolant cvadratic neted precum cel folosit în regula Simpson va da duce la rezultate bune. Oricum, este des întâlnit cazul ca funcţia pe care încercăm să o integrăm este netedă peste interval. Aceasta înseamnă că orice funcţia este puternic oscilantă, ori ea nu are derivate în anumite puncte. În aceste cazuri, regula Simpson poate da rezultate foarte slabe. O cale uzuală de a rezolva această problemă este spargerea intervalului într-un număr de subintervale mici. Apoi, regula Simpson este aplicată pe fiecare subinterval, se sumează rezultatele pentru o obţine o aproximare a integralei pe întreg intervalul. Acest tip de abordare se numeşte regula Simpson compusă. Presupunem că intervalul este împărţit în n subintervale, cu n un număr par. Atunci, regula Simpson compusă este dată de:

97

,

unde pentru cu ; în particular, şi

. Formula de mai sus poate fi scrisă astfel:

Eroarea făcută de regula Simpson compusă este mărginită (în valoarea absolută) de

, unde .

Această formulare împarte intervalul în subintervale de lungime egală. În practică, este deseori avantajos de a folosi subintervale de lungimi diferite, şi de a ne concentra eforturile în locurile unde integrantul se comportă mai puţin bine. Aceasta conduce la metoda Simpson ajustată.

98

7.2.6. Metoda ajustativă Simpson Metoda ajustativă Simpson, numită de asemenea regula ajustativă Simpson, este o metodă de integrare numerică propusă de William M. McKeeman în 1962. Este, probabil, primul algoritm ajustativ pentru integrarea numerică care a apărut, deşi mai multe metode moderne bazate pe integrarea Gauss – Kronrod şi integrarea Clenshaw – Curtis sunt, în general, preferate acum. Metoda ajustativă Simpson foloseşte o estimare a erorii pe care o obţinem calculând o integrală definită folosind regula lui Simpson. Dacă eroare depăşeşte o toleranţă specificată de utilizator, algoritmul cere subdivizarea intervalului de integrare în două şi aplicarea metodei ajustative a lui Simpson pe fiecare subinterval într-o manieră recursivă. Tehnica este de obicei mult mai eficientă decât regula de compunere Simpson din moment de foloseşte mai puţine funcţii de evaluare în punctele în care funcţia este bine aproximată de o funcţie integrală. Criteriul de determinare a momentului în care oprim subdivizarea intervalului este:

unde este un interval care are mijlocul c, şi sunt estimările date de regula lui Simpson pe intervalele corespunzătoare şi este toleranţa dorită pentru interval. Regula Simpson este un caz particular al metodei Romberg. Teoria acestei metode arată că regula lui Simpson este exactă atunci când integrantul este un polinom de grad 3 sau mai mic. Potrivit metodei Romberg, cea mai corectă estimare Simpson pentru şase valori ale funcţiei se combină cu cea mai puţin corectă estimare Simpson

pentru trei valori ale funcţiei aplicând corecţia

.

Aici, constanta 15 este aleasă pentru a asigura că după ce aplicăm corecţia, este obţinută o estimare care este exactă pentru polinoame de grad cel mult 5.

7.3. Integrarea în puncte neechidistante Pentru funcţii care nu sunt periodice, oricât, de departe, cele mai precise metode cu puncte neechidistante sunt cvadratura gaussiană şi cvadratura Clenshaw – Curtis.

7.3.1. Cvadratura gaussiană În analiza numerică, o regulă de integrare este o aproximare a integralei definite a unei funcţii, uneori formulată ca suma ponderată a

99

valorilor funcţiei în punctele date înăuntrul domeniului de integrare. O regula de cvadratură gaussiană n – puncte, numită după Carl Friedrich Gauss, este in regulă de integrare construită pentru a produce un rezultat exact pentru polinoamele de grad , printr-o alegere corespunzătoare a n punctelor xi

şi n ponderilor wi. Domeniul de integrare pentru o astfel de regulă este, prin convenţie, luat [−1, 1], astfel că regula este exprimată ca

.

Se poate arăta (vezi Press, ş.a. sau Stoer şi Bulirsch) că punctele de evaluare sunt chiar rădăcinile polinomului care aparţine clasei polinoamelor ortogonale. Reguli pentru problema de bază Pentru problema de integrare exprimată mai sus, polinoamele asociate sunt polinoamele Legendre, . Cu cel de-al n – lea polinom normalizat pentru a obţine Pn(1) = 1, nodul Gauss i, xi, este a i – a rădăcină a lui Pn; ponderea ei este dată de (Abramowitz & Stegun 1972, p. 887)

.

În continuare sunt date câteva reguli de integrare pentru ordine mici:

Numărul de puncte, n

Punctele, xi Ponderile, wi

1 0 2

2 1

3

0

4

100

5

0

Schimbarea intervalului pentru cvadratura gaussiană O integrală pe [a, b] trebuie schimbată într-o integrală pe [−1, 1] înainte să aplicăm regula de integrare gaussiană. Această schimbare a intervalului poate fi făcută astfel:

.

După aplicarea regulii de integrare Gauss se obţine următoarea aproximare:

Alte forme ale cvadraturii gaussiene Problema integrării poate fi exprimată într-o manieră puţin mai generală prin introducerea unei funcţii pondere , pozitivă, în integrant, şi permiţând un interval diferit de [−1, 1]. Adică, problema este calculul

pentru unele alegeri ale lui a, b şi ω. Pentru a = −1, b = 1 şi ω(x) = 1, problema este aceeaşi cu cea considerată mai sus. Alte alegeri conduc la alte reguli de integrare, unele dintre ele fiind cele din tabelul următor.

Intervalu

lω(x) Polinoamele ortogonale

[−1,1] 1 Polinoamele Legendre

(−1,1) Polinoamele Jacobi

101

(−1,1)Polinoamele Chebyshev (de prima

speţă)

[−1,1]Polinoamele Chebyshev (de speţa a

doua)

[0,∞) Polinoamele Laguerre

(−∞,∞) Polinoamele Hermite

Teorema fundamentală

Fie q un polinom netrivial de grad n astfel încât

pentru toţi .Dacă alegem nodurile să fie zero-urile lui q, atunci există ponderile wi

care fac integrala calculată exactă pentru toate polinoamele de grad 2n − 1 sau mai mic. În plus, toate aceste noduri se vor găsi în intervalul deschis (a, b).

Estimarea eroriiEroarea unei reguli de integrare Gauss poate fi exprimată după cum

urmează. Pentru un integrant care are 2n derivate continue,

pentru un ξ oarecare din (a, b), unde pn este polinomul ortogonal de grad n şi

unde

În cazul special ω(x) = 1 avem estimarea erorii

.

Stoer şi Bulirsch au remarcat că nu este convenabil în practică, din moment ce poate fi dificil de estimat derivata de ordinul 2n, şi, în plus, eroarea actuală poate fi mult mai mică decât o margine stabilită de derivată. O altă abordare este de a folosi regulile de cvadratură gaussiană de ordine diferite, şi să estimăm eroarea ca diferenţa dintre cele două rezultate. Pentru acest scop, pot fi utile regulile de integrare Gauss – Kronrod.

102

Regulile Gauss – Kronrod Dacă intervalul [a,b] este împărţit, punctele de evaluare Gauss ale noilor subintervale coincid cu punctele de evaluarea anterioare (exceptând zero, pentru numerele impare), şi astfel, integrantul trebuie să fie evaluat în fiecare punct. Regulile Gauss – Kronrod sunt extinderi ale regulilor de cvadratură gaussiană generate prin adăugarea a n + 1 puncte la o regulă n – punct în aşa fel încât regula rezultantă să fie de ordin 3n + 1. Acest fapt permite calculul estimărilor de ordin superior refolosind valorile funcţiei a unei estimări de ordin mic. Diferenţa dintre regula de cvadratură gaussiană şi extensia ei Kronrod este des utilizată ca o estimare a erorii aproximării. Regulile sunt numite după Alexander Kronrod care le-a inventat în 1960. Algoritmii în QUADPACK au la bază regulile Gauss–Kronrod. Un exemplu popular combină o regulă Gauss 7 – puncte cu o regulă Kronrod 15 – puncte. Deoarece punctele Gauss sunt încorporate în punctele Kronrod, un total de doar 15 evaluări produc estimarea integrării cât şi eroarea estimării.

Nodurile Gauss Ponderile

±0.94910 79123 42759 0.12948 49661 68870

±0.74153 11855 99394 0.27970 53914 89277

±0.40584 51513 77397 0.38183 00505 05119

0.00000 00000 00000 0.41795 91836 73469

Nodurile Kronrod Ponderile

±0.99145 53711 20813 0.02293 53220 10529

±0.94910 79123 42759 0.06309 20926 29979

±0.86486 44233 59769 0.10479 00103 22250

±0.74153 11855 99394 0.14065 32597 15525

±0.58608 72354 67691 0.16900 47266 39267

±0.40584 51513 77397 0.19035 05780 64785

±0.20778 49550 07898 0.20443 29400 75298

0.00000 00000 00000 0.20948 21410 84728

Patterson a arătat cum se găsesc alte extensii de acest tip.7.3.2. Cvadratura tanh – sinh Cvadratura tanh – sinh este o metodă de integrare numerică introdusă de Hidetosi Takahasi şi Masatake Mori în 1974. Ea foloseşte schimbarea de variabilă

103

pentru a transforma o integrală pe într-o integrală pe întreaga dreaptă reală. Pentru un pas dat h, integrala este aproximată de suma

cu: - abscisele

şi - ponderile:

.

La fel ca cvadratura gaussiană, cvadratura tanh-sinh este potrivită pentru integrarea cu precizie arbitrară, unde este dorită o acurateţe de sute sau mii de digiţi.

Convergenţa este pătratică pentru integranţi cu comportare suficient de bună: dublarea numărului de puncte de evaluare dublează grosolan numărul de digiţi corecţi. Cvadratura tanh – sinh este mai puţin eficientă decât cvadratura Gauss pentru integranţi netezi, dar, spre deosebire de cvadratura Gauss lucrează bine şi cu integranţi care au singularităţi sau derivate infinite într-unul sau în ambele capete ale intervalului de integrare. Un alt avantaj este acela că abscisele şi ponderile sunt relativ simplu de calculat. Costul calculului perechilor abscise – ponderi pentru o acurateţe de n – digiţi este

comparabil cu pentru cvadratura Gauss. Comparând schema cu cvadratura Gauss şi cvadratura funcţiei eroare, Bailey ş.a. au găsit că schema tanh – sinh „pare a fi cea mai bună pentru integranţi de tipul celor mai des întâlniţi în cercetările matematice experimentale”.

7.3.3. Cvadratura Clenshaw – Curtis Cvadratura Clenshaw-Curtis şi cvadratura Fejér sunt metode de integrare numerică, sau de cvadratura numerică, bazate pe o dezvoltare a integrantului în termeni de polinoame Chebyshev. Echivalent, ele folosesc o schimbare de variabilă x = cosθ şi folosesc o aproximare a transformării cosinus discretă pentru seriile cosinus. În afară că au o convergenţă rapidă

104

comparabilă cu regulile de integrare gaussiană, cvadratura Clenshaw – Curtis şi Fejér conduc natural către regulile de integrare (în care diferite ordine de precizie împart punctele), care este importantă pentru ambele integrări (multidimensionale şi ajustativă). Pe scurt, funcţia care trebuie integrată este evaluată în punctul de extrem N sau în rădăcinile unui polinom Chebyshev şi aceste valori sunt folosite pentru a construi o aproximare polinomială pentru această funcţie. Această aproximare este atunci integrată exact. În practică, ponderile de integrare pentru valoarea funcţiei în fiecare nod sunt precalculate, şi acest calcul poate fi executat în timpul O(NlogN) cu ajutorul diferitor algoritmi legaţi de transformării rapide Fourier pentru DCT. O cale simplă de a înţelege algoritmul este de a înţelege faptul că cvadratura Clenshaw – Curtis (propus de acei autori în 1960) se referă la o integrare printr-o schimbare de variabilă x = cosθ. Algoritmul este exprimat pentru integrala unei funcţii pe un interval (orice alt interval poate fi obţinut printr-o rescalare). Pentru această integrală, putem scrie

.

Aceasta înseamnă că trebuie să transformăm această problemă din integrarea lui în integrarea lui . Acest lucru poate fi făcut dacă ştim seria cosinus pentru :

caz în care integrala devine:

.

Bineînţeles, pentru a calcula coeficienţii seriei cosinus

105

trebuie să executăm încă o dată o integrare numerică, deci pentru început aceasta poate să pară că nu simplifică problema. Spre deosebire de calculul integralelor arbitrare, oricum, integrările seriilor Fourier pentru funcţii periodice (precum , prin construcţie), până la frecvenţa Nyquist k =

N, sunt precis calculate de punctele N echidistante şi ponderile

echidistante pentru (cu excepţia punctelor terminale care sunt

ponderate de , pentru a evita dubla numărare). Adică, aproximăm integrala

seriilor cosinus prin transformarea cosinus discretă (DCT) de tipul I:

pentru şi apoi folosim formula de mai sus pentru integrală în

termenii acelor . Legătura cu polinoamele Chebyshev Motivul pentru care aceasta este legată de polinoamele Chebyshev

este acela că, prin definiţie, şi astfel, seria cosinus

de mai sus este într-adevăr o aproximare a lui prin polinoame Chebyshev:

,

şi astfel integrăm “cu adevărat” integrând dezvoltarea aproximantă în

termenii polinoamelor Chebyshev. Punctul de evaluare

corespunde extremei polinomului Chebyshev . Faptul că o astfel de aproximare Chebyshev este chiar o serie cosinus după o schimbare de variabile este responsabilă pentru convergenţa rapidă a aproximării deoarece mai mulţi termeni sunt incluşi. O serie cosinus converge foarte repede pentru funcţiile pare, periodice, şi suficient de netede. Acest lucru este adevărat aici, din moment de este pară şi periodică în prin construcţie, şi este de k ori derivabilă în orice punct dacă este de k ori derivabilă pe . (Spre deosebire, aplicând direct o dezvoltare în serie cosinus lui în loc de în mod uzual nu va converge rapid deoarece panta dezvoltărilor pare – periodice va fi, în general, discontinuă).

106

Comparaţia cu cvadratura Gauss Metoda clasică a integrării Gauss evaluează integrantul în puncte şi este construită pentru a integra exact polinoamele de grad maxim

. În contrast, cvadratura Clenshaw – Curtis, mai întâi evaluează integrantul în puncte şi integrează exact polinoamele doar până la gradul N. Poate părea, deci, că cvadratura Clenshaw – Curtis este intrinsec, mai rea, decât cvadratura Gauss, dar, în realitate acesta nu pare să fie cauza. În practică, mai mulţi autori au observat că cvadratura Clenshaw – Curtis poate avea o precizie comparabilă cu cea a cvadraturii gaussiene pentru acelaşi număr de puncte. Acest lucru este posibil deoarece cei mai mulţi integranţi numerici nu sunt polinoame (în special, de când polinoamele pot fi integrate analitic), şi aproximarea mai multor funcţii în termeni ai polinoamelor Chebyshev converge rapid. De fapt, rezultate teoretice recente (Trefethen, 2006) argumentează că atât cvadratura gaussiană cât şi cea Clenshaw – Curtis au eroarea mărginită de

pentru un integrant de k ori diferenţiabil.Un avantaj des citat al integrării Clenshaw – Curtis este acela că

ponderile cvadraturii pot fi evaluate în timpul prin algoritmi rapizi de transformare Fourier (şi analogii lor pentru DCT), pe când

ponderile integrării gaussiene cer un timp de calcul.

Ca o chestiune practică, o cvadratura numerică de ordin ridicat este de puţine ori făcută prin simpla evaluare a formulei de integrare pentru N foarte mare.

În schimb, de obicei foloseşte o schemă de integrare ajustativă care mai întâi evaluează integrala de ordin mic, şi apoi, succesiv, rafinează precizia doar în regiunile unde integrala este inexactă. Pentru a evalua precizia integralei, comparăm răspunsul cu cel al regulii de integrare pentru ordine mici pare. Ideal, această regulă de integrare de ordin mic evaluează integrantul pe o submulţime a punctelor N originale, să minimizez evaluările integranţilor. Aceasta se numeşte regulă de integrare nested, şi aici cvadratura Clenshaw – Curtis are avantajul că regula pentru ordinul N foloseşte o submulţime a punctelor de la ordinul 2N. În contrast, regulile de integrare gaussiene nu sunt natural nested, şi astfel trebuie să folosească cvadratura Gauss – Kronrod sau metodele similare.

107

7.3.4. Cvadratura Fejér Fejér a propus două reguli de integrare foarte asemănătoare cu cele de integrare Clenshaw – Curtis, dar cu mult timp înainte (în 1933). A “doua” regulă de integrare a lui Fejér este identică cu cea dată de Clenshaw – Curtis. Singura diferenţă este aceea că punctele terminale şi sunt fixate în zero. Adică, Fejér a folosit doar extremele interioare ale polinoamelor Chebyshev, adică punctele staţionare. “Prima” regulă de integrare a lui Fejér evaluează estimând

într-o altă mulţime de puncte echidistante, la jumătatea distanţei

dintre extreme: pentru . Acestea sunt rădăcinile

lui , cunoscute sub denumirea de nodurile Chebyshev. (Aceste mijloace echidistante sunt singurele celelalte alegeri ale punctelor de integrare care păstrează atât simetria pară a transformării cosinus cât şi simetria translaţională a seriei Fourier periodice) . Aceste lucruri conduc la:

care este chiar DCT de tip II. În ciuda faptului că Fejér a descoperit aceste tehnici înainte lui Clenshaw şi Curtis, denumirea “cvadratura Clenshaw – Curtis” a devenit cea standard.

7.3.5. Cvadratura ajustativă În matematica aplicată, cvadratura ajustativă este un proces în care integrala unei funcţii este aproximată folosind regulile de integrare statice pe subintervale adaptativ rafinate ale domeniului de integrare. În general, algoritmii ajustativă sunt la fel de eficienţi şi eficaci precum cei tradiţionali pentru integranţi “cu comportare bună”, dar, sunt de asemenea eficaci pentru integranţi “cu comportare rea”, pentru care algoritmii tradiţionali eşuează. Schema generală

Cvadratura ajustativă urmează schema generală:1. procedura integrate (f,a,b,tau)

2.

3.

4. dacă atunci

108

5. m = (a + b) / 26. Q = integrate(f,a,m,tau) + integrate(f,m,b,tau)7. sfârşit dacă8. return Q Este calculată o aproximare Q a integralei lui pe intervalul (linia 2), precum şi o eroare estimată (linia 3). Dacă eroarea estimată este mai mare decât toleranţa cerută (linia 4), atunci intervalul este împărţit (linia 5) şi integrala este aplicată pe ambele intervale separat (linia 6). Se returnează atât estimarea iniţială sau suma înjumătăţită calculată recursiv (linia 7). Componentele importante sunt însăşi regula de integrare

,

estimarea erorii

,

şi logica deciderii împărţirii cărui interval, şi apoi când să ne oprim. Bineînţeles că există mai multe variante ale acestei scheme. Cea mai uzuală va fi discutată mai târziu. Reguli de integrare de bază Regulile de integrare au în general forma:

unde nodurile şi ponderile sunt, în general, calculate anterior. În cazul general, sunt folosite formulele Newton – Cotes de ordin par, unde nodurile sunt uniform spaţiate în interval:

.

Când sunt folosite asemenea reguli, punctele în care f(x) a fost evaluată pot fi reutilizate după revenire:

109

O strategie similară este folosită la cvadratura Clenshaw - Curtis, unde nodurile sunt alese astfel:

.

sau, când este folosită cvadratura Fejér:

.

Un algoritm poate fi ales să folosească diferite metode de integrare pe subintervale diferite, de exemplu, folosind o metodă de ordin înalt doar acolo unde integrantul este neted. Estimarea erorii Unii algoritmi de integrare generează o secvenţă de rezultate care ar trebui să se apropie de valoarea corectă. Altfel, putem folosi o “regulă nulă” care are forma regulii de integrare de mai sus, dar ale cărei valori vor fi zero pentru un integrant simplu (de exemplu, dacă integrantul este un polinom de grad apropiat). Subdiviziunea logică Cvadratura ajustativă “locală” face eroarea acceptabilă pentru un interval dat proporţional cu lungimea acelui interval. Acest criteriu poate fi cu greu satisfăcut dacă integrantul se comportă rău în doar câteva puncte, de exemplu. Pentru câţiva paşi de discontinuitate. Alternativ, putem cere doar ca suma erorilor pe fiecare subinterval să fie mai mică decât cererea utilizatorului. Aceasta va fi cvadratura ajustativă “globală”. Cvadratura ajustativă globală poate fi mai eficientă (folosind mai puţine evaluări ale integrantului) dar este, în general, mai dificil de programat şi poate cere mai mult spaţiu de lucru pentru înregistrarea informaţiei pe mulţimea curentă de intervale.

1x 2x 3x

1x 2x 3x

2x 3x

adân

cim

ea

1x

110

7.4. Integrări cu funcţii pondere Mai general, putem pune problema integrării unei funcţii în locul unei funcţii pondere care este cunoscută de dinainte:

.

Cazul cel mai comun este , ca mai sus, dar în unele aplicaţii este utilă o funcţie pondere diferită. Motivul fundamental este că, din moment ce poate fi luat în socoteală aprioric, eroarea de integrare poate fi făcută să depindă doar de precizia în aproximarea lui , indiferent de cât de rău se comportă funcţia pondere. Cvadratura Clenshaw – Curtis poate fi generalizată la acest caz după cum urmează. Ca şi înainte, ea merge găsind dezvoltarea cos-sin a funcţiei via DCT, şi apoi integrând fiecare termen a seriei cosinus. Acum, aceste integrale sunt de forma:

.

Pentru majoritatea funcţiilor , această integrală nu poate fi calculată analitic, spre deosebire de cazul anterior. Din moment ce aceeaşi funcţie pondere este folosită în general pentru mai mulţi integranţi ,

putem permite să calculăm aceşti numeric pentru a mări precizia în prealabil. În plus, din moment ce w(x) este specificată în mod general analitic, putem uneori folosi metode speciale să calculăm Wk. De exemplu, au fost dezvoltate diferite metode pentru a aplica cvadratura Clenshaw – Curtis pentru integranţi de forma f(x)w(x) cu o funcţie pondere w(x) puternic oscilantă, de exemplu o sinusoidă sau o funcţie Bessel. Acest lucru este util pentru precizia ridicată a calculului seriilor Fourier sau a seriilor Fourier – Bessel, în care metodele de integrare pentru w(x) = 1 sunt problematice din cauza preciziei ridicate cerută să rezolve contribuţia oscilaţiilor rapide. Aici, partea integrantului cu oscilaţii rapide este luată în seamă prin metode speciale pentru Wk, pe când funcţia necunoscută f(x) se comportă mai bine. Alt caz în care funcţiile pondere sunt utile, în special, dacă integrantul este necunoscut dar are o singularitate cunoscută de o anumită formă, de exemplu, o discontinuitate cunoscută sau o divergenţă a integralei (precum

) într-un anumit punct.

111

De notat este că cvadratura gaussiană poate fi de asemenea adaptată pentru diferite funcţii pondere, dar tehnica este întrucâtva diferită. În cvadratura Clenshaw – Curtis, integrantul este întotdeauna evaluat în aceeaşi mulţime de puncte a lui w(x), corespunzătoare extremelor sau rădăcinilor polinomului Chebyshev. În cvadratura gaussiană, diferite funcţii pondere conduc la diferite polinoame ortogonale, şi astfel, diferite rădăcini în care integrantul este evaluat.

7.5. Metoda Nyström Metoda Nyström de discretizare a unei ecuaţii integrale foloseşte o regulă de integrare; adică, aplicând regula de integrare

de exemplu, ecuaţiei Fredholm neomogene de speţa a doua

rezultă

.

7.6. Formula Euler-MacLaurin În matematică, formula Euler – MacLaurin furnizează o legătură puternică între integrale şi sume. Ea poate fi folosită pentru aproximarea integralelor prin sume finite, sau invers, pentru evaluarea sumelor finite şi seriilor infinite folosind integralele şi calculul. Formula a fost descoperită independent de Leonhard Euler şi Colin MacLaurin în jurul anului 1735 (şi a fost generalizată mai târziu ca formula lui Darboux). Euler a avut nevoie de ea pentru a calcula seriile infinite încet convergente în timp ce MacLaurin a folosit-o pentru a calcula integrale. Formula Dacă n este un număr natural şi este o funcţie netedă (în sensul, diferenţiabilă) definită pentru toate numerele reale x dintre 0 şi n, atunci

integrala poate fi aproximată de suma:

112

(vezi regula trapezului). Formula Euler – MacLaurin furnizează expresii

pentru diferenţa dintre sumă şi integrală în termeni ai derivatelor în 0 şi

n. Pentru orice număr natural p, avem:

unde B1 = −1/2, B2 = 1/6, B3 = 0, B4 = −1/30, B5 = 0, B6 = 1/42, B7 = 0, B8 = −1/30, ... sunt numerele Bernoulli, şi R este eroarea care este mică pentru valori corespunzătoare ale lui p. Utilizând substituţia, putem adapta această formulă de asemenea la funcţii f care sunt definite pe alte intervale ale dreptei reale. Restul Restul R este dat de

,

unde sunt polinoamele Bernoulli. Restul poate fi estimat prin:

.

Sume care implică un polinom Dacă f este un polinom şi p este suficient de mare, atunci restul dispare. De exemplu, dacă f(x) = x3, putem alege p = 2 pentru a obţine, după simplificare,

.

Integrarea numerică Formula Euler – MacLaurin este, de asemenea, folosită pentru analiza detaliată a erorilor în integrarea numerică; în particular, metodele de extrapolare depind de ea. Dezvoltări asimptotice ale sumelor În contextul calculului dezvoltărilor asimptotice ale sumelor şi seriilor, de obicei, cea mai utilă formă a formulei Euler – MacLaurin este:

,

unde a şi b sunt întregi. Adesea, dezvoltarea rămâne validă chiar şi după ce trecem la limită pentru sau , sau amândouă. În multe cazuri integrala sau membru drept pot fi evaluate în forma închisă în termenii funcţiilor elementare chiar şi atunci când suma din membrul stâng nu poate.

113

Atunci toţi termenii seriilor asimptotice pot si exprimaţi în termenii funcţiilor elementare. De exemplu,

.

Aici, membrul stâng este egal cu suma dintre şi , unde,

aceasta din urmă este funcţia poligamma de prim ordin definită prin

; funcţia gamma este egală cu dacă z este

un întreg pozitiv. Scăzând din ambele părţi rezultă o dezvoltare

asimptotică pentru . Dezvoltarea, în schimb, serveşte ca punct de

plecare pentru una din derivaţiile estimării precise a erorii pentru aproximarea Stirling a funcţiei factorial. Derivaţia din inducţie matematică Polinoamele Bernoulli Bn(x), n = 0, 1, 2, ... pot fi definite recursiv după cum urmează:

,

şi pentru .

Primii termeni sunt:

,

Valorile Bn(1) sunt numerele Bernoulli. Pentru n ≥ 2, avem Bn(0) = Bn(1). Funcţiile Bernoulli periodice Pn sunt date de

pentru ,

adică, coincid cu polinoamele Bernoulli pe (0,1) şi sunt periodice, de perioadă 1. Fie integrala

,

114

unde

(din moment ce ),

. Integrând prin părţi, obţinem

Sumăm pentru k de la 1 la şi rezultă

.

Adunând în ambii membri şi rearanjând, avem

(1)

Ultimii doi termeni dau eroarea atunci când integrala este folosită pentru aproximarea sumei. Considerăm

,

unde

Integrând din nou, prin părţi, obţinem,

115

Sumând apoi pentru k de la 1 la , şi apoi înlocuind ultima integrală în (1) cu ce am arătat că este egală, avem:

.

Procedând în continuare la fel se obţine demonstraţia prin inducţie matematică a formulei de sumare Euler – MacLaurin, în care, pasul de inducţie se bazează pe integrarea prin părţi şi pe identităţile pentru funcţiile periodice Bernoulli. Pentru o obţine limitele mărimii erorii atunci când suma este aproximată de integrală, notăm că polinoamele Bernoulli pe intervalul [0,1] îşi ating valorile maxime absolute în extreme şi Bn(1) este al n – lea număr Bernoulli. Derivaţia din analiza funcţională Formula Euler – MacLaurin poate fi înţeleasă ca o curioasă aplicaţie a unor idei din spaţiile Hilbert şi analiza funcţională. Fie Bn(x) polinoamele Bernoulli. Funcţiile duale polinoamelor Bernoulli sunt date de

unde δ este funcţia delta a lui Dirac. Deasupra este o notaţie formală pentru ideea de a lua derivatele într-un punct; astfel avem:

pentru n > 0 şi o funcţie arbitrară, dar diferenţiabilă, f(x) definită pe intervalul unitate. Pentru cazul n =0, definim . Polinoamele Bernoulli, împreună cu dualele lor, formează o mulţime ortogonală a intervalului unitate: avem

şi

.

Formula de sumare Euler – MacLaurin rezultă atunci ca o integrală peste cea din urmă. Avem

116

Luând x = 0, şi rearanjând termenii, obţinem formula cunoscută, împreună cu eroarea. De notat că numerele Bernoulli sunt definite ca Bn = Bn(0), şi că acestea se anulează pentru n impar mai mare decât 1. Această derivaţie nu presupune că f(x) este netedă şi că se comportă bine; anume, că f poate fi aproximată prin polinoame; echivalent, că f este o funcţie analitică reală. Formula de sumare Euler – MacLaurin poate fi astfel văzută a fi un rezultat al reprezentării funcţiei pe intervalul unitate prin produsul direct al polinoamelor Bernoulli şi al dualelor lor. Oricât, reprezentarea nu este completă pe mulţimea funcţiile pătrat integrabile. Dezvoltarea în termenii polinoamelor Bernoulli are un nucleu netrivial. În particular, sin(2πnx) este într-un nucleu; integrala lui sin(2πnx) se anulează pe intervalul unitate, fiind diferenţa derivatelor ei în capetele intervalului.

7.7. T – integrarea T – integrarea este o tehnică de integrare numerică dezvoltată de Jon Michael Smith în 1970 pentru a facilita comanda şi controlul avionului. Prescurtare pentru “integrare numerică tunabilă” ea foloseşte un pas fix ca mărime şi o formulă iterativă care depinde de faza şi parametrii câştigaţi. Fie integrantul şi P şi G faza şi parametrii câştigaţi. În plus,

limita din membrul stâng al integralei este notată prin şi este mărimea pasului. T – integrarea se defineşte prin formula recursivă:

.

Aici, este notaţia pentru . Cantitatea aproximează

.

Dacă G = 1, atunci metoda se reduce la următoarele tehnici binecunoscute de integrare numerică pentru valorile date de P:

P = 0: regula stângă a dreptunghiului, P = 1/2: regula trapezului, P = 1: regula dreaptă a dreptunghiului.

117

T – integrarea poate fi acordată la problema care este folosită pentru rezolvare. T – integrarea este bazată pe teoria informaţiei, nu pe teoria aproximării. T – integrarea are un domeniu simplu de frecvenţă de ajustare a parametrilor: un parametru de ajustare a fazei şi un parametru de ajustare câştigat. Interesant, pentru probleme open-loop, fixând câştigul şi variind faza obţinem toţi integratorii numerici de prim ordin şi un infinit de noi integratori până acum necunoscuţi. Pentru aplicaţiile closed loop T – integrarea conduce la o infinitate de integratori neclasici care produce integrarea numerică exactă a sistemelor liniare şi integrarea aproape exactă a sistemelor neliniare. Pentru aplicaţiile asupra sistemelor de informaţie (calculator, control, comunicare şi simulare) T – integratorul simplu de prim ordin face toţi integratorii numerici bazaţi pe teoria clasică a aproximării. Simularea mişcării avionului pentru diferite configuraţii aviatice (roată sus, roată jos, clapă sus, clapă jos, motor afară etc.) şi condiţiile dinamice (decolarea, aterizarea etc.) devine o simplă problemă de tunare a T – integratorului la condiţiile de zbor care sunt simulate. În acest sens, T – integratorul se adaptează la problema pe care încearcă să o rezolve.

BIBLIOGRAFIE

[1] Bucur,C.M., Metode numerice, Ed. Facla, Timişoara, 1973.[2] Ciurea,E., Algoritmi Introducere în algoritmica grafurilor, Ed. Tehnică, Bucureşti, 2001.[3] Coman,G., Analiză numerică, Ed. Libris, Cluj, 1995.[4] Croitoru,C., Tehnici de bază în optimizarea combinatorie, Ed. Universităţii „Al. I. Cuza”, Iaşi, 1992.[5] Cuculescu,I., Analiză numerică, Ed. Tehnică, Bucureşti, 1967.[6] Demidovici,B.P., Maron,I., Elements de calcul numerique,

118

Ed. Mir de Mosou, 1973.[7] Dodescu,Gh., Toma,M., Metode de calcul numeric, E. D. P., Bucureşti, 1976.[8] Dodescu,Gh., Metode numerice în algebră, Ed. tehnică, Bucureşti, 1979.[9] Ichim,I., Marinescu,G., Metode de aproximare numerică, Ed. Academiei R. S. R., Bucureşti, 1986.[10] Ignat,C., Ilioi,C., Jucan,T., Elemente de informatică şi calcul numeric, Univ. „Al. I. Cuza”, Iaşi, Fac. de Matematică, 1989.[11] Juan Antonio Infante del Rio, Jose Maria Rey Cabezas, Metodos Numericas,Teoria,problemas y practicas con

MATLAB, Ed. Piramide, 2002.[12] Knut,D.E., Sortare şi căutare, vol. 3, Tratat de programarea calculatoarelor, Ed. Tehnică, Bucureşti, 1976.

[13] Livovschi,L., Georgescu,H., Sinteza şi analiza algoritmilor, Ed. Ştiinţifică şi Enciclopedică, Bucureşti, 1986.[14] Melhorn,K., Data Structures and Algorithms, Springer-Verlag, Berlin, 1984.[15] Mihu,C., Metode numerice în algebra liniară, Ed. Tehnică, Bucureşti, 1977.[16] Popovici,P., Cira,O., Rezolvarea numerică a ecuaţiilor neliniare, Ed. Signata, Timişoara, 1992.[17] Press,W.H., Teuklosky,S.A., Vetterling,W.T., Flannery,B.P., Numerical Recipes in C: The Art of scientific Computing, (Cambridge University Press, Cambridge, 1992).[18] Scheiber,E., Metode numerice, Univ. Transilvania din Braşov, Facultatea de Matematică – Informatică, (electronic).[19] Toma,M., Odăgescu,I.,

119

Metode numerice şi subrutine, Ed. Tehnică, Bucureşti, 1980.[20] Vladislav,T., Raşa,I., Analiză numerică, Ed. Tehnică, Bucureşti, 1997.[21] Vraciu,G., Popa,A., Metode numerice cu aplicaţii în tehnica de calcul, Scrisul românesc, Craiova, 1982.[22] ***, Borland International, Inc. Borland C++, Programming Guide (Borland International, Scotts wally, CA, 1992).

120