Analiza numerica.pdf curs de matematica

7/29/2019 Analiza numerica.pdf curs de matematica

http://slidepdf.com/reader/full/analiza-numericapdf-curs-de-matematica 1/96

UNIVERSITATEA ”OVIDIUS” CONSTANTA

FACULTATEA DE MATEMATICA SI INFORMATICA

Prof. univ. dr. CONSTANTIN POPA Lect. univ. dr. ELENA PELICAN

ANALIZA NUMERICA

Constanta 2006



Introducere

Lucrarea de fata ısi propune sa prezinte principalele tehnici de aproxi-

mare din analiza numerica clasica. Ea are la baza cursurile, seminariilesi laboratoarele din acest domeniu tinute de autori ın perioada 1991-2004la specializarile de matematica, informatica si inginerie din UniversitateaOvidius din Constanta. Bagajul de cunostinte matematice necesare par-curgerii lucrarii se rezuma la cursurile de analiza matematica , algebraliniara si ecuat ii diferentiale ordinare care se tin, ın mod uzual ın anii I siII (semestrul I) la toate specializarile mentionate mai sus. Prezentarea re-spectivilor algoritmi este ınsotita de un minimum de consideratii teoreticecare asigura atat ıntelegerea proprietatilor si comportarii lor ın aplicatii catsi o baza pentru abordari ulterioare ale unor extinderi si generalizari. Inlucrare au fost prezentate metodele clasice din diverse capitole ale analizeinumerice, insistandu-se pe sublinierea ideilor fundamentale legate de algo-ritmii respectivi. Astfel, abordarea unor variante moderne si ımbunatatiteva fi mai usor de realizat pentru cititorii interesat i. Seturile de exercitii careıncheie fiecare capitol contin complemente, dezvoltari si aplicatii practiceale metodelor descrise. Detalii legate de anumite demonstratii, extinderiale metodelor prezentate sau solutii ale unor exercitii se pot gasi ın lucrarile[22], [23], [24], [21], [20].

i



Notatii

IN multimea numerelor naturaleZ multimea numerelor ıntregiIR multimea numerelor realeC multimea numerelor complexeIRn, Cn spatiile euclidiene n dimensionale real,

respectiv complexIN ∗, Z∗, IR∗, C∗ IN ∗ = IN \ 0; analog pentru celelalteIR[X ] multimea polinoamelor ın nedeterminata X,

cu coeficienti numere realegradP gradul polinomului P C n([a, b]) multimea functiilor definite pe [a, b], de n ori derivabile,

cu derivata de ordin n continuaC ∞([a, b]) multimea functiilor definite pe [a, b], indefinit derivabileMn(K) multimea matricelor patrate de dimensiune n ≥ 1

cu elemente din corpul K,unde K este IR sau C; pentru K = IR vom notaprin Mn

ai sau (A)i linia i din matricea A ∈ Mn

aij sau (A)ij elementul de pe pozitia (i, j) din matricea A ∈ Mn.

zi sau (z)i componenta de pe pozitia i din vectorul z ∈ IRn

detA determinantul matricei A ∈ Mn

tr(A) urma matricei A ∈ Mn (tr(A) = a11 + · · · + ann)At transpusa matricei A

∈ Mn

diagd1, . . . , dn matrice diagonala cu elementele d1, . . . , dn pediagonala principala

diag(A) matricea diagonala ce are pe diagonala principalaelementele de pe diagonala matricei A

k(A) numarul de conditionare al matricei inversabile A ınraport cu o norma data · ;k(A) = A · A−1

σ (A) spectrul unei matrice A ∈ Mn (valorile proprii)ρ (A) raza spectrala a unei matrice A ∈ Mn

(ρ(A) = max|λ|, λ ∈ σ(A))δ ij simbolul lui Kronecker

O(hn) marime de ordin hn



Cuprins

1 Reprezentarea numerelor ın calculator. Erori de calcul 11.1 Reprezentarea p-adica a numerelor reale . . . . . . . . . . . 11.2 Reprezentarea numerelor reale ın calculator. Erori de rotun-

jire si trunchiere . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Exercitii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Calculul valorilor functiilor elementare. 122.1 Schema lui Horner . . . . . . . . . . . . . . . . . . . . . . . 122.2 Schema lui Horner generalizata . . . . . . . . . . . . . . . . 142.3 Dezvoltari ın serie Taylor . . . . . . . . . . . . . . . . . . . 16

3 Aproximarea solutiilor ecuatiilor si sistemelor de ecuatiineliniare 213.1 Metoda bisectiei . . . . . . . . . . . . . . . . . . . . . . . . 21

3.2 Metoda coardei . . . . . . . . . . . . . . . . . . . . . . . . . 233.3 Metoda aproximatiilor succesive pe IR . . . . . . . . . . . . 263.4 Metoda lui Newton pe IR . . . . . . . . . . . . . . . . . . . 283.5 Metoda aproximatiilor succesive si metoda Newton pentru

sisteme de ecuatii neliniare . . . . . . . . . . . . . . . . . . 303.6 Exercitii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4 Aproximare, interpolare si derivare numerica 354.1 Aproximarea cu polinoame Bernstein . . . . . . . . . . . . . 364.2 Interpolare cu polinoame Lagrange . . . . . . . . . . . . . . 374.3 Functii spline cubice . . . . . . . . . . . . . . . . . . . . . . 394.4 Derivarea numerica . . . . . . . . . . . . . . . . . . . . . . . 404.5 Exercitii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

5 Integrare numerica 455.1 Formule de cuadratura de tip Newton-Cotes . . . . . . . . . 46

5.1.1 Aproximarea integralelor pe IR . . . . . . . . . . . . 465.1.2 Aproximarea integralelor pe IR2 . . . . . . . . . . . . 53

ii



5.2 Formule de cuadratura cu noduri Gauss . . . . . . . . . . . 555.2.1 Formule de tip Cebısev . . . . . . . . . . . . . . . . 55

5.2.2 Formule bazate pe polinoame Legendre . . . . . . . 575.2.3 Formule de tip Gauss pentru integrale duble . . . . . 57

5.3 Exercitii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

6 Aproximarea solutiilor ecuatiilor diferentiale ordinare 626.1 Metode de tip Taylor . . . . . . . . . . . . . . . . . . . . . . 636.2 Metode de tip Runge-Kutta . . . . . . . . . . . . . . . . . . 666.3 Exercitii . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

7 Metode directe pentru sisteme liniare 707.1 Metoda eliminarii a lui Gauss . . . . . . . . . . . . . . . . . 707.2 Descompunerea LU . . . . . . . . . . . . . . . . . . . . . . . 74

7.2.1 Conditii suficiente pentru ca EG sa functioneze farapivotare . . . . . . . . . . . . . . . . . . . . . . . . . 77

7.3 Descompunerea Choleski . . . . . . . . . . . . . . . . . . . . 797.4 Exercitii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

8 Aproximarea valorilor si vectorilor proprii 828.1 Metoda Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . 838.2 Metoda puterilor . . . . . . . . . . . . . . . . . . . . . . . . 858.3 Metoda puterilor inverse . . . . . . . . . . . . . . . . . . . . 878.4 Exercitii . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Bibliografie 89



Capitolul 1

Reprezentarea numerelor ın

calculator. Erori de calcul

1.1 Reprezentarea p-adica a numerelor reale

Definitia 1 Fie (an)n ⊂ IR si S n =n

k=1

ak, n ≥ 1. Dac˘ a sirul (S n)n≥1

este convergent la un num˘ ar x ∈ IR spunem c˘ a seria n≥1

an este convergent˘ a

si are suma x si scriem n≥1

an = x. (1.1)

Atunci S n se numeste suma part ial˘ a de rang n, iar an termenul general (de rang n) al seriei.

Definitia 2 Dac˘ a n≥1

an = x este o serie convergent˘ a definim cantitatea

Rn = x − S n, n ≥ 1 (1.2)

numit˘ a restul de ordin n al seriei. Evident c˘ a

limn→∞Rn = 0. (1.3)

Propozitia 1 (i) Dac˘ a n≥1

an este convergent˘ a, atunci limn→∞an = 0.

(ii) Avem relat ia

Rn = limk→∞

n+k j=n+1

a j

. (1.4)

1



(iii) Dac˘ a

n≥1

an si

n≥1

bn sunt serii de numere reale pozitive, seria

n≥1

bn

este convergent˘ a si exist a n0 ≥

1 si c > 0 astfel ınc at

an ≤ c · bn, ∀ n ≥ n0, (1.5)

atunci si seria n≥1

an e convergent˘ a.

Corolarul 1 Fie p ∈ IN, p ≥ 2 si (an)n ⊂0, 1, ..., p-1. Atunci seria n≥1

an

pneste convergent˘ a.

Observatia 1 Orice num˘ ar real a se poate exprima ın scrierea zecimal˘ a

a = a0, a1a2a3...... (1.6)

ın care a1, a2,... se numesc cifre zecimale ( ai ∈ 0, 1, . . . , 9, ∀ i ∈ IN ),iar a0 este partea ıntreag˘ a a lui a (notat˘ a ın continuare prin [a]). Chiar si pentru a ∈ Q ın scrierea (1.6) pot ap˘ area o infinitate de cifre zecimale (numere periodice). Scrierea ”pozit ionala” (1.6) este echivalent˘ a din punct de vedere matematic cu

a = a0 +n≥1

an

10n, (1.7)

seria din dreapta fiind convergent˘ a ın baza corolarului anterior.

Observatia 2 Scrierea (1.6) nu este unic˘ a ∀ a ∈ IR. De exemplu pentru a = 1/2 avem

a = 0.5 = 0.49999.... (1.8)

Rezultatul care urmeaza generalizeaza consideratiile anterioare pentru oricebaza p ∈ IN, p ≥ 2 (pentru demonstratii vezi [22]).

Teorema 1 (i) Orice num˘ ar real α ∈ [0, ∞) admite o reprezentare de forma

α = a0 +n≥1

an

pn, (1.9)

cu an ∈ 0, 1,..,p − 1, ∀ n ≥ 1.(ii) Scrierea (1.9) este unic˘ a dac˘ a si numai dac˘ a

α = k

pn, ∀ k, n ∈ IN. (1.10)

2



Corolarul 2 Toate afirmat iile din enunt ul Teoremei 1 r˘ amˆ an valabile si pentru numere α

∈(

−∞, 0).

Observatia 3 Din Teorema 1 rezult˘ a c˘ a scrierea (1.9) nu afecteaz˘ a va-loarea (ca num˘ ar real) a p˘ art ii ıntregi sau p˘ art ii zecimale a num˘ arului α.Reprezentarea (1.9) se va nota prin a = a0, a1a2a3......( p) (vezi (1.6)).

Pentru numere ıntregi obtinem o scriere similara cu (1.9), data de urmatorulrezultat.

Corolarul 3 Fie α ∈ Z. Atunci exist˘ a n ≥ 1 (care depinde de α) si b0, b1, ..., bn ∈ 0, 1,...,p − 1 unic determinat i astfel ınc at

α = ±n

k=0

bk · pk. (1.11)

Definitia 3 Fie p ∈ N, p ≥ 2. Definim

Q p = pk · n | k, n ∈ Z, (1.12)

numit˘ a multimea numerelor practice (ın baza p).

Teorema 2 Mult imea Q p este dens˘ a ın IR.

1.2 Reprezentarea numerelor reale ın calculator.

Erori de rotunjire si trunchiere

Vom prezenta pentru ınceput cateva chestiuni cu privire la reprezenta-rea numerelor reale ın virgula mobila (vezi pentru detalii [19], [25]). Pentrunumere ıntregi, a ∈ Z, reprezentarea se face exact, ın baza p = 2, pe unanumit numar n de biti ce depinde de sistemul de calcul folosit. Vom exem-plifica ın continuare pentru cazul n = 8. Astfel, pentru a = 5, reprezentareaar fi (conform corolarului 3)

0 0 0 0 0 1 0 1

Daca a este negativ, de exemplu a = −5, se utilizeaza metoda coduluicomplementar, adica ın locul lui a se reprezinta numarul

a = 28 + a. (1.13)

In cazul nostru a = 28 − 5 = 2 5 6 − 5 = 251 va avea reprezentarea1 1 1 1 1 0 1 1

3



Acest tip de reprezentare, desi exacta, impune limitari ale valorilor nu-merelor ce pot fi ınregistrate. Intr-adevar avem

| a | ≤ 2n−1 − 1. (1.14)

In acest fel toate numerele ıntregi pozitive, respectiv negative vor avea0, respectiv 1 pe prima pozitie din stanga (ceea ce permite implicit re-cunoasterea semnului lor). Fie acum a ∈ IR \ Z, a > 0. Din Teorema 1(pentru p = 2) rezulta

a = a0 +n≥1

an

2n. (1.15)

Daca a0 > 0, fie m cel mai mare numar natural cu proprietatea ca

2m

≤a

0(1.16)

si n0 dat den0 = m + 1. (1.17)

Atunci, din (1.11), (1.16), (1.15) si (1.17) rezulta

a = (n≥1

an2n

) · 2n0. (1.18)

De asemenea, din definitia lui n0 obtinem ca

a1= 0. (1.19)

Daca a0 = 0, atunci definim direct n0 ca fiind cel mai mare numar negativcu proprietatea ca ın scrierea (1.18) este ındeplinita conditia (1.19).

Figura 1.1: Reprezentarea pe 32 biti

Exemplul 1 Pentru a = 0, 000101001(2) avem n0 = −3 si a = 2−3 ·0, 101001(2). Cantit˘ at ile

M a =n≥1

an2n

(1.20)

4



si C a = 64 + n0 (1.21)

se numesc mantisa (normalizata), respectiv caracteristica num˘ arului a. Dac˘ a a ∈ IR \ Z, a < 0 se utilizeaz˘ a schema anterioar˘ a pentru b =−a > 0. Astfel, reprezentarea unui num˘ ar a ∈ IR \ Z (pe 32 de bit i) se va face ca in Figura 1.1.

Observatia 4 Reprezentarea caracteristicii pe 7 bit i impune restrict ia

0 ≤ C a ≤ 127 (1.22)

sau

−64

≤n0

≤63 (1.23)

ceea ce determin˘ a limit˘ ari ale valorilor maxime si minime ale lui a.

Observatia 5 Deoarece num˘ arul de bit i afectat i mantisei este finit rezult˘ a c˘ a numerele cu care lucreaz˘ a un calculator electronic constituie o submult ime a mult imii numerelor practice din sect iunea 1.1.

Limitarea (fizica) a numarului de pozitii pentru mantisa determina aparitiaerorilor de trunchiere si de rotunjire (la introducerea datelor reale sau ınurma operatiilor aritmetice elementare). Regulile uzuale de transformare aunui numar real ce se utilizeaza de obicei sunt urmatoarele: daca numarula are ın mantisa mai mult de 24 de cifre binare ai, i ≥ 1 atunci

1. daca a25

= 0 se ınregistreaza a1

, a2

,...,a24

neschimbate.

2. daca a25 = 1 si ak = 0, ∀ k ≥ 26 se procedeaza ca ın cazul 1.

3. daca a25 = 1 si ∃ k ≥ 26 astfel ıncat ak = 1 se adauga o unitate laa24 si se opereaza modificarile respective ın mantisa (spre stanga).

Erorile ce apar ın cazurile 1. si 2. se numesc erori de trunchiere, iar celedin cazul 3. erori de rotunjire. Aparent, ele se ”concentreaza” asuprapozitiei 24 din mantisa, dar valoarea absoluta a numarului a poate fi afec-tata chiar ın pozitii din stanga virgulei dupa cum vom vedea ın exempleleurmatoare.

Exemplul 2 Fie a = 1001, 11 . . . 1 21 pozit ii

ın baza 2. ˆ In scrierea (1.18) vom avea

n0 = 4, a25 = 1 si ak = 0, ∀ k ≥ 26. Atunci, conform regulii 2 de mai sus,a va fi trunchiat la valoarea

a = 1001, 11 . . . 1 20 pozit ii

5



Exemplul 3 Fie a = 1001, 11 . . . 1

22 pozit ii

ın baza 2. ˆ In scrierea (1.18) vom avea

n0 = 4, a25 = 1, a26 = 1, ak = 0, ∀ k ≥ 27.

Conform regulii 3, a va fi rotunjit la valoarea a = 1010, 00....0.

Exemplul 4 Fie a = 1 . . . 1 24 pozit ii

0100, 01 ın baza 2. ˆ In scrierea (1.18) vom

avea n0 = 28, a25 = 1. Conform regulii 1, a se trunchiaz˘ a la valoarea

a = 1 . . . 1 24 pozit ii

0000, 00....

Exemplul 5 Fie a = 1 . . . 1 25 pozit ii

0 1 0 0, 0 1 ın baza 2. ˆ In scrierea (1.18)

vom avea n0 = 29, a25 = 1, a26 = 0 dar a27 = 1. Atunci, conform regulii 3,a se va rotunji la valoarea

a = 1 0 . . . 0 24 pozit ii

00000, 00....

Pentru a ”masura” efectul acestor rotunjiri si trunchieri vom introducecateva notiuni specifice.

Definitia 4 Fie a, a ∈ IR. Vom numi eroare absoluta de aproximare a lui a si vom nota cu a cantitatea

∆(a) =| a − a | . (1.24)

Dac˘ a a = 0, cantitatea

δ (a) =∆(a)

| a | (1.25)

se va numi eroare relativa de aproximare a lui a prin a. Vom numi cifre semnificative toate cifrele din scrierea zecimal˘ a (ın baza 10) a unui num˘ ar, ıncep and cu prima diferit˘ a de zero din stˆ anga. O cifr˘ a semnificativ˘ a a aproxim˘ arii a a lui a se va numi exacta dac˘ a eroarea absolut˘ a ∆(a) nu depaseste unitatea ordinului respectivei cifre.

Exemplul 6 Fie numerele

a = 2, 003450; a = 2, 0034400.

Atunci ∆(a) = 0, 00001,

toate cifrele lui a sunt semnificative, dar exacte sunt doar 2, 0, 0, 3 si primul 4.

6



Observatia 6 Pentru exemplele 2 − 5 de mai ınainte, erorile absolute de aproximare sunt urm˘ atoarele(ın ordine).

Exemplul 2 ∆(a) = 2−21

Exemplul 3 ∆(a) = 2−22

Exemplul 4 ∆(a) = 22

Exemplul 5 ∆(a) = 23 + 2−2

Valorile corespunzatoare exemplelor 4 si 5 sunt incomparabil mai mari(chiar ıngrijoratoare!) decat cele din exemplele 2 si 3. Totusi, trebuie satinem cont si de valoarea absoluta a numerelor din exemplele 4 si 5. Intr-adevar, daca am calcula ın aceste cazuri erorile relative am obtine valoricam de acelasi ordin cu cele absolute corespunzatoare exemplelor 2 si 3.Acest fapt explica si necesitatea introducerii notiunii de eroare relativa. Eaeste un fel de ”eroare procentuala” care arata ce parte reprezinta eroareaabsoluta din ıntregul numar (adica eroarea absoluta ∆(a) = 8, 25 din ex-emplul 5 este mare teoretic vorbind, dar ea nu afecteaza esential numarul acare are o valoare absoluta foarte mare ≈ 1010). Acelasi lucru se ıntamplasi la numere foarte mici, dar ın sens invers.

Exemplul 7 Fie num˘ arul

a = 0, 0 . . . 0 24 pozit ii

1 0 . . . 0 23 pozit ii

101

ın baza 2. ˆ In scrierea (1.18) vom avea n0 =−

24, a25

= 1, a27

= 1, deci conform regulii 3, a se va rotunji la valoarea

a = 0, 0 . . . 0 24 pozit ii

1 0 . . . 0 22 pozit ii

1.

ˆ In acest caz avem ∆(a) = 2−51 · 3 si δ (a) = 2−25

Exista diverse variante de studiu teoretic al propagarii erorilor de calcul(vezi [25], [19] si referintele lor). Acestea ınsa nu fac obiectul lucrarii defata. Trebuie totusi sa punctam cateva aspecte importante:

a) Toate aceste studii prezinta ın final majorari pentru erorile absolute,respectiv relative ce apar ın urma introducerii datelor si a operatiilor ele-mentare. In plus, cel putin din punct de vedere teoretic, se demonstreazaca aceste majorari sunt optimale, ın sensul ca exista situatii cand ele suntatinse. De exemplu, ın cazul sumei a doua numere a, b ∈ IR, daca a, b suntdoua aproximari ale lor, se poate demonstra ca (vezi [25]).

∆(a + b) ≤ ∆(a) + ∆(b) (1.26)

7



si ca exista a, b astfel ıncat ın (1.26) sa avem egalitate.b) In urma acestor studii teoretice apar si unele reguli de ”aranjare” a

calcului astfel ıncat cumularea erorilor de calcul sa fie minima. De exempluıntr-o suma este mai bine sa grupam mai ıntai termenii dupa ordinul demarime si/sau dupa semn. Acestea ınsa au ın gene-ral un efect relativ. Esential este ca algoritmii care se utilizeaza sa fieoptimizati din punct de vedere al eficientei (vezi sectiunile urmatoare).

c) Retinem totusi ideea propagarii si acumularii de erori de calcul ceeace ındeamna la o anumita doza de prudenta ın abordarea etapei finale aunei probleme de analiza numerica, aceea a calcului efectiv al aproximarilorrespective.

d) Spre deosebire de multimea numerelor reale care este infinita si fara”goluri” /”spatii libere” ıntre numere (adica ıntre oricare doua numere reale

diferite, exista ıntotdeuna un altul si pentru orice numar real dat exista unaltul mai mare si unul mai mic decat el) multimea numerelor reprezentate ınvirgula mobila (”machine numbers”) este finita (existand un cel mai mare/cel mai mic numar (pozitiv)) si este discreta. Astfel, multimea numerelorreale normalizate ın virgula mobila se defineste ca

F B,t = M · BE /M,E ∈ Z , M = 0 sau Bt−1 ≤ |M | < BtB - se numeste baza si este un numar natural mai mare sau egal cu 2, t -numarul de zecimale, M - mantisa, iar E - exponentul.

Multimea numerelor reprezentate ın calculator (pe scurt, numere masina)se defineste astfel:

MB,t,α,β = g ∈ F B,t, α ≤ E ≤ β unde B , t α , β - sunt date de constructia sistemului de calcul (nu suntmemorate explicit si nu pot fi schimbate). Pentru fiecare numar de masinase memoreaza doar valorile M si E care reprezinta de fapt numarul M ·BE .

In formatul IEEE (Institute for Electrical and Electronic Engineers)adoptat din 1985 de majoritatea producatorilor de microcomputere, sefoloseste B = 2, iar precizia poate fi single (t = 24), double (t = 53)sau extended (t = 64) (de exemplu, pentru reprezentarea numerelor ex-tended, primul bit este de semn, urmeaza 11 biti pentru exponent/caracteristicasi ultimii 52 de biti pentru mantisa).

Doua numere masina consecutive sunt g = M · BE si g = (M + 1)BE ,iar distanta relativa dintre ele este:

g − g

g=

1

M ≤ 1

Bt−1

In MB,t,α,β exista:

8



• cel mai mare numar masina pozitiv, L = (Bt − 1)Bβ

• cel mai mic numar masina pozitiv, s = B

t

−1

· B

α

Toate numerele care apar ın calcule si sunt ın intervalul (−s, s) nu pot fireprezentate si drept urmare, sunt setate cu 0 (mesa j de ”underflow”), iarnumerele mai mari ca L sau mai mici ca −L determina oprirea efectuariicalculelor (mesaj de ”overflow”).

In formatul IEEE,t α β s L

Single 24 −149 104 1.18E − 038 3.40E + 038Double 53 −1074 971 2.23E − 308 21024

Extended 64 −16445 16320 2−16382 216384

e) In conformitate cu standardele IEEE, rezultatul unei operat ii arit-

metice elementare a (+, −, ·, /) este cel mai ”corect” posibil ın sensul ca este”rotunjit” la cel mai apropiat numar masina. Daca ”·” si ”/” nu sunt criticedin punct de vedere al erorii, ”+” si ”−” pot produce surprize din cauzaanularii bitilor semnificativi datorate scaderii numerelor foarte apropiate cavaloare. De exemplu, pentru a = 2.145648xx si b = 2.145611xx si folosindreprezentarea pe 6 biti, valoarea exacta pentru c = a − b este 0.000037xx;valoarea lui a ın reprezentarea pe 6 biti este a = 0.214564xxx iar a lui beste b = 0.214561xxx si deci c = a −b = 0.000003xxx, ceea ce conduce laeroarea relativa

|c −

c|

|c

|=

|0.000037 − 0.000003|

|0.000037

|=

0.000034

0.000037=

34

37= 0.(918).

Urmatorul exemplu se refera la situatii ın care asociativitatea adunarii numai are loc. Astfel, ın reprezentarea pe 6 biti, avem

(1 + 1000000) − 1000000 = 1000000 − 1000000 = 0

iar1 + (1000000 − 1000000) = 1 − 0 = 1.

1.3 Exercitii

1. Fie integrala I n = 1

0 xnex

10dx pentru n natural. Pentru n = 0,

valoarea integralei este usor de calculat si este egala cu I 0 = 10e0.1 −10 1.05170918075647624811707826490.

Pentru n > 0, integrala se poate calcula folosind urmatoarea formularecursiva (rezultata din integrarea prin parti)

I n = 10e0.1 − 10nI n−1.

9



Din pacate ınsa, folosind aceasta formula, obtinem rezultate care nusunt conforme cu realitatea, acest algoritm fiind numeric instabil.

(a) Scrieti un program C care calculeaza I n pentru n = 0, 1,.., 6si foloseste rotunjirea la p = 3, 6, 9, 12 zecimale ın calcule. Deasemenea programul trebuie sa afiseze o coloana cu va-lori cat de exacte posibil, adica calculate fara extra rotunjiri.Cazul p = 0 semnifica faptul ca nu s-a aplicat nici o extra rotunjire (vezi functia de rotunjire de mai jos, pentru detalii). Tabelulva trebui sa aiba formatul de mai jos (cu valorile corespunzatoaredate de programul de calcul).

p=0 p=3 p=6 p=9 p=12

n=0 1.051709e+00 1.100000e+00 1.051700e+00 1.051709e+00 1.051709e+00

n=1 5.346174e-01 1.00000e-01 5.347000e-01 5.346172e-01 5.346174e-01n=2

n=3

n=4

n=5

n=6

Chiar daca rotunjirea cu o anumita precizie nu este realizata deo functie standard a limbajului de programare ales, ea poate fiobtinuta, de exemplu, cu ajutorul unei rutine de tipul urmator.

double eround(double x, int p)

int d;double temp;

if (x == 0|| p <= 0) return x;

/*nu face nimic daca p ≤ 0(nu rotunjeste) sau x = 0*/

/*(ar calcula lg 0, nu convine)*/

d = (int)(log 10(fabs(x)));

if (fabs(x) > 1) p − −;

/*ajusteaza p ( p se poate modifica ın functie de */

/*variabila care este apelata prin valoarea ei */

/* absoluta)*/

temp = rint(x ∗ ttpo( p − d));

/*modificare si rotunjire x */return ((double)temp ∗ ttpo(d − p));

/*remodificare si returnare valoare*/

Aceasta rutina are urmatorul efect: eround(x,p) ıntoarce val-oarea lui x rotunjita la p zecimale seminificative. Va trebui saapelati aceasta functie de fiecare data cand o operatie ın virgula

10



mobila ar fi putut sa schimbe numarul zecimalelor seminifica-tive, avand ca efect rotunjirea noii valori la p ze-

cimale semnificative.(b) Calculati cat mai exact urmatoarea integrala:

J n =

1

0xne

−x

10000dx

pentru n = 6. Explicati tehnica folosita si rezultatele obtinute.

2. Fie y = x − sin x. Se stie ca pentru valori mici ale lui x, sin x xsi deci scaderea lor poate duce la ”surprize” neplacute. Cum poate fiacest lucru evitat?Indicat ie. Se dezvolta ın serie Taylor sin x (vezi sectiunea 2.3) si se

folosesc formulelet1 =

x3

6

tn+1 = − tnx2

(2n + 2)(2n + 3), n ≤ 1.

Astfel, notand sn =n

k=1

tk, obtinem

s1 = t1

sn+1 = sn + tn+1.

3. Aproximati e−5 prin:

(a) e−5 9

i=0

(−5)i

i!=

9

i=0

(−1)i · 5i

i!

(b) e−5 =1

e5 1

9i=0

5!

i!

Valoarea aproximativa de referinta pentru e−5 ın acest exercitiu este6 · 74 × 10−3. Care din formulele de mai sus da un rezultat mai bunsi de ce?

11



Capitolul 2

Calculul valorilor functiilor

elementare.

2.1 Schema lui Horner

Fie P (x) = a0xn + a1xn−1 + ... + an, n ≥ 1 un polinom cu coeficientireali si ξ ∈ R fixat. Nu este indicat ca pentru valori mai mari ale lui n,calculul valorii P (ξ ) sa se efectueze direct ın expresia lui P (x) ( datoritaerorilor de rotunjire ce pot aparea la calculul puterilor xn). Modalitateauzuala de calcul se bazeaza pe un algoritm simplu, dar eficient, care ınlaturaneplacerile mentionate anterior - schema lui Horner data de urmatoareasecventa de operatii aritmetice elementare

b0 = a0, bi = ai + bi−1ξ, 1 ≤ i ≤ n, P (ξ ) = bn. (2.1)

Observatia 7 Pe lˆ ang˘ a simplitatea schemei de calcul remarc˘ am si faptul c˘ a b0,...,bn−1 reprezint˘ a coeficient ii cˆ atului ımp art irii lui P la x − ξ (ın ordinea descresc˘ atoare a puterilor lui x) bn = P (ξ ) fiind restul ımp˘ art irii.

De multe ori ın algoritmii de calcul este nevoie ca, dupa o schimbare devariabila de forma x = y + ξ sa se determine coeficientii noului polinom ınvariabila y.

Q(y) = P (y + ξ ) = A0yn + ... + An (2.2)

Algoritmul utilizat ın acest sens este urmatorul: Pentru y = 0, din relatia

P (x) = Q(y) = P (y + ξ ) (2.3)

rezulta

An = P (ξ ). (2.4)

12



Fie P 1(x) polinomul dat de

P (x) = (x − ξ )P 1(x) + P (ξ ) (2.5)

Pentru y = x − ξ , din (2.2) rezulta

P (x) = (x − ξ )[A0(x − ξ )n−1 + ... + An−1] + P (ξ ) (2.6)

de unde, folosind unicitatea lui P 1, obtinem

P 1(x) = A0(x − ξ )n−1 + ... + An−1 (2.7)

si, pentru x = ξ An−1 = P 1(ξ ). (2.8)

Rationamentul se reia punand

P 1(x) = (x − ξ )P 2(x) + P 1(ξ ). (2.9)

Astfel, coeficientii An, An−1,...,A0 din (2.2) vor fi dati de relatiile

Ak = P n−k(ξ ), k = n, n − 1,..., 0 (2.10)

unde am notat P 0(x) = P (x), iar polinoamele P 1,...,P n sunt cele obtinutesuccesiv ın schema anterioara.

Exemplul 8 Fie P (x) = x4 − 8x3 + 5x2 + 2x − 7 si ξ = 2. Algoritmul prezentat anterior poate fi scris concentrat ca in Figura 2.1. Astfel, ın

raport cu (2.2) si (2.10) obt inem

Q(y) = P (y + 2) = y4 − 19y2 − 42y − 31.

Figura 2.1: Schema lui Horner

13



2.2 Schema lui Horner generalizata

Conform teoremei de ımpartire cu rest, dat fiind polinomul Q cu coefici-enti reali si grad(Q) ≥ 1 exista doua polinoame cu coeficienti reali C ,respectiv R, cu grad(R) < grad(Q) astfel ıncat

P (x) = C (x)Q(x) + R(x). (2.11)

In cazul grad(Q) = 1, conform observatiei din sectiunea 2.1, putem deter-mina coeficientii lui C si valoarea lui R (ın acest caz R = 0 sau grad(R) = 0)fara a efectua ımpartirea, cu ajutorul schemei lui Horner. Aceeasi problemase poate pune ınsa si pentru grad(Q) ≥ 2. Avem doua situatii distincte.

Cazul 1. Q are toate radacinile reale. Atunci problema se poaterezolva prin aplicari succesive ale algoritmului din sectiunea anterioara.

Vom ilustra acest lucru printr-un exemplu. Fie Q(x) = (x − a)(x − b). Din(2.11) rezulta

P (x) = (x − a)(x − b)C (x) + (αx + β ) (2.12)

si pentru a determina coeficientii lui C , α si β aplicam schema lui Hornerrelativ la ımpartirea lui P la x − a. Obtinem

P (x) = (x − a)C 1(x) + β 1 (2.13)

Repetam apoi pasul anterior cu C 1 ın locul lui P si x − b ın locul lui x − a,

obtinand C 1(x) = (x − b)C 2(x) + α1, (2.14)

de unde, folosind (2.13) rezulta

P (x) = (x − a)(x − b)C 2(x) + α1(x − a) + β 1. (2.15)

Dar deoarece polinoamele C si R din (2.11) sunt unic determinate, din(2.15) obtinem

C (x) = C 2(x), α = α1, β = β 1 − aα1. (2.16)

Extinderea la cazul general

Q(x) = (x − a1)...(x − an) (2.17)

este imediata.

Cazul 2. Q are si radacini complexe. Deoarece Q are coeficientireali aceste radacini sunt perechi de forma a ± ib. Problema revine atunci

14



la a construi un algoritm care sa permita determinarea coeficientilor luiC si α, β din relatia

P (x) = (x2 + b1x + b0)C (x) + αx + β (2.18)

fara a efectua ımpartirea propriu-zisa (ın (2.18) factorul x2 + b1x + b0 areradacinile a ± ib). Vom construi ın continuare un asemenea algoritm, nu-mit schema lui Horner generalizata, plecand de la relatiile ce apar ıntrecoeficientii polinoamelor din (2.18). Pentru simplitatea expunerii, vom con-sidera cazul particular P (x) = a3x3 + a2x2 + a1x + a0. In urma efectuariiımpartirii, obtinem ın (2.18)

C (x) = a3x + (a2 − a3b1), (2.19)

α = a1 − a3b0 − b1(a2 − a3b1), (2.20)

β = a0 − b0(a2 − a3b1). (2.21)

Notandb1 = −b1, b0 = −b0, (2.22)

formam urmatorul tablou

a3 a2 a1 a0

0 0 a3b0 b0(a2 + a3)0 a3b1 b1(a2 + a3b1) 0

(2.23)

Observam ca suma termenilor de pe coloanele 1, 2, 3, 4 reprezinta chiarcoeficientii lui C (x) si α, β din (2.18)-(2.21). O analiza mai atenta neconduce la urmatorul algoritm de constructie a coeficientilor din tabloul(2.23). Sa consideram mai ıntai tabloul initial

a3 a2 a1 a0

0 0 α13 α14

0 α22 α23 0(2.24)

ın care elementele α22, α13, α23, α14 sunt necunoscute. Ele se determina ınfelul urmator:α

22=suma elementelor de pe coloana 1 ınmultita cu b

1;

α13=suma elementelor de pe coloana 1 ınmultita cu b0;α23=suma elem. de pe coloana 2 ınmultita cu b1;α14 = suma elem. de pe col.2 ınmultita cu b0.Deducem astfel urmatorul algoritm - dat fiind polinomul

P (x) = anxn + an−1xn−1 + ... + a0,

15



se efectueaza urmatorii pasi:Pasul 1. Se formeaza tabloul urmator

an an−1 . . . an− j+1 . . . a1 a0

0 0 . . . α1 j . . . α1n α1,n+1

0 α22 . . . α2 j . . . α2n 0

Σ1 = an Σ2 . . . Σ j . . . Σn Σn+1

(2.25)

ın care cunoastem an, an−1, . . . , a1, a0, iar Σ2, . . . , Σn+1 reprezinta sumeleelementelor de pe coloanele 2, . . . , n + 1 din tabel.Pasul 2. Coeficientii αij se calculeaza astfel

α22 = an · (−b1),

iar pentru j = 3, . . . , n α1 j = Σ j−2 · (−b0)α2 j = Σ j−1 · (−b1)

(2.26)

α1,n+1 = Σn−1(−b0).

Pasul 3. Valorile an, Σ2, . . . , Σn−1 reprezinta coeficientii puterilor lui x (ınordine descrescatoare) ale catului C (x) din (2.18), iar α, β sunt date de

α = Σn, β = Σn+1 (2.27)

Observatia 8 Combinˆ and cele dou˘ a scheme de tip Horner (clasic˘ a si ge-

neralizat˘ a) putem determina cˆ atul si restul ımpart irii lui P la orice polinom unitar Q. ˆ In plus, ambii algoritmi sunt usor de programat.

2.3 Dezvoltari ın serie Taylor

Definitia 5 Fie I = (α, β ) ⊂ IR, f : I −→ IR o funct ie de clas˘ a C ∞, n ∈N∗ si ξ ∈ I fixat. Expresia

T n(x; ξ ) = f (ξ ) +f (ξ )

1!(x − ξ ) + . . . +

f (n)(ξ )

n!(x − ξ )n, x ∈ I (2.28)

se numeste polinomul Taylor de grad n ın jurul lui ξ asociat funct iei f ,

iar Rn(x; ξ ) = f (x) − T n(x, ξ ) (2.29)

restul Taylor de ordin n. De asemenea vom nota cu n≥0

f (n)(ξ )

n!(x − ξ )n (2.30)

16



sau

f (ξ ) +f (ξ )

1!

(x

−ξ ) + . . . +

f (n)(ξ )

n!

(x

−ξ )n + . . . (2.31)

seria Taylor ın jurul lui ξ a funct iei f .

Se cunosc urmatoarele rezultate (vezi pentru detalii [26]).

Teorema 3 (Formula Taylor cu restul Lagrange)Pentru f , I , ξ ca mai ınainte, pentru orice x ∈ I , exist˘ a δ > 0 cuprins ıntre x si ξ astfel ınc at

f (x) = T n(x; ξ ) +(x − ξ )n+1

(n + 1)!f (n+1)(δ ). (2.32)

Corolarul 4 (Formula McLaurin)ˆ In ipotezele teoremei 1, dac˘ a 0 ∈ I , atunci

f (x) = T n(x; 0) +xn+1

(n + 1)!f (n+1)(δ ), x ∈ I (2.33)

Teorema 4 (Reprezentarea ın serie Taylor)Fie f : I = (α, β ) −→ IR de clas˘ a C ∞ si [a, b] ⊆ I presupunem c˘ a ∃ M > 0astfel ınc at

f (n)(x)

≤ M, ∀ x ∈ [a, b], ∀ n ≥ 0. (2.34)

Atunci,

∀ξ

∈(a, b) seria Taylor a lui f ın jurul lui ξ converge uniform pe

[a, b] la f . Vom nota aceasta prin

f (x) =n≥0

f (n)(ξ )

n!(x − ξ )n, x ∈ [a, b] (2.35)

sau

f (x) = f (ξ ) +f (ξ )

1!(x − ξ ) + . . . +

f (n)(ξ )

n!(x − ξ )n + . . . (2.36)

Observatia 9 Din (2.29) si (2.32) obt inem pentru Rn(x; ξ ) exprimarea

Rn(x; ξ ) =(x

−ξ )n+1

(n + 1)! f (n+1)(δ ), (2.37)

cu x ∈ I si δ ıntre x si ξ .

Observatia 10 ˆ In (2.32) ın locul lui δ se poate pune ξ + θ(x − ξ ),cu θ ∈ (0, 1).

17



Observatia 11 Funct iile elementare admit exprim˘ ari de tipul (2.32) sau (2.33) sau dezvolt˘ ari de tipul (2.36) pe anumite intervale I

⊂IR. Aceste

dezvolt˘ ari pot fi utilizate pentru calculul aproximativ al valorilor acestor funct ii. Vom prezenta ın continuare aceste lucruri ın cazul cˆ atorva funct ii elementare uzuale.

Functia exponentialaFunctia exponentiala f : IR −→ (0, ∞), f (x) = ex admite dezvoltarea

(de tip (2.33))

ex = 1 + x + . . . +xn

n!+ Rn(x), (2.38)

valabila pentru orice x ∈ IR cu Rn(x) dat de

Rn =

eθx

(n + 1)! , θ ∈ (0, 1). (2.39)

Pentru x ∈ [−1, 1] din (2.38) si (2.39) rezulta

|ex − T n(x; 0)| = |Rn(x)| ≤ e

(n + 1)!, (2.40)

care reprezinta o buna evaluare a priori a erorii de aproximare a valoriiexacte ex prin valoarea polinomului Taylor T n(x; 0). Pentru |x| ≥ 1 avem

x = [x] + q, q ∈ [0, 1) (2.41)

si

ex

= e[x]

· eq

. (2.42)Pentru e[x] procedam astfel: folosind (2.38), cu evaluarea (2.40) obtinem obuna aproximare a lui e sau 1/e (dupa cum x > 0 sau x < 0) si cu aceastacalculam e[x] prin ınmultiri succesive. Pentru eq se utilizeaza direct (2.38).

Observatia 12 ˆ In evaluarea valorilor polinomului Taylor din (2.38), pen-tru a evita calculul factorialelor calcul˘ am T n(q ; 0) prin

T n(q ; 0) = u0 + u1 + . . . + un, (2.43)

unde ui se obt in recursiv prin

u0 = 1, uk =q

·uk−

1

k , k = 1, . . . , n (2.44)

Functia logaritmicaPentru x ∈ (−1, 1] avem

ln(1 + x) = x − x2

2+

x3

3− x4

4+ . . . + (−1)n+1 xn

n+ . . . (2.45)

18



De asemenea, pe acelasi interval, ln(1 + x) admite si o reprezentare de tipul(2.28)

−(2.29) ın care restul de ordin n este dat de

Rn(x) =xn+1

n + 1· (−1)n+1

(1 + θx)n+1, θ ∈ (0, 1). (2.46)

Exista ınsa urmatoarele dezavantaje ın ceea ce priveste utilizarea dezvoltarii(2.45) cu restul dat de (2.46) pentru aproximarea valorilor functiei ln(1+x) :

1. domeniul limitat de valori x ∈ (−1, 1] ⇒ 1 + x ∈ (0, 2];

2. pentru |x| ”apropiat” ca valoare de 1 se stie ca seria (2.45) convergefoarte ”ıncet”, iar evaluarea erorii data de (2.46) conduce la valori

mult prea mari ale lui n (pentru o buna aproximare).

Aceste inconveniente pot fi eliminate prin cateva transformari efectuateasupra dezvoltarii (2.45) (care se bazeaza pe convergenta uniforma a serieirespective). Astfel, pentru x ∈ (−1, 1) din (2.45) rezulta

ln(1 − x) = −x − x2

2− x3

3− . . . − xn

n− . . . (2.47)

Prin scaderea relatiilor (2.47) si (2.45) obtinem

ln(1 − x

1 + x) =

−2(x +

x3

3+

x5

5+ . . .) (2.48)

Pentru z = 1−x1+x , cum x ∈ (−1, 1), obtinem z ∈ (0, ∞) si din (2.48)

ln z = −21 − z

1 + z

+

1

3

1 − z

1 + z

3+ . . .

. (2.49)

Observatia 13 Dezvoltarea (2.49) este valabil˘ a pentru orice z ∈ (0, ∞),deci tot domeniul de definit ie al funct iei ln z.

Vom analiza acum restul ın dezvoltarea (2.49). Stim ca ∀ z > 0,∃ m ∈ IN, t ∈ [ 1

2 , 1], unic determinati, cu proprietatea

z = 2m · t. (2.50)

Atunci

ln z = ln(2m · t) = m · ln 2 + ln t = m · ln 2 − 2n

k=1

ξ 2k−1

2k − 1− Rn(ξ ), (2.51)

19



cu ξ ∈ (0, 1/3] dat de

ξ =1 − t

1 + t

. (2.52)

Observand ca2

1 − ξ 2≤ 9

4, (2.53)

rezulta pentru Rn(ξ ) urmatoarea evaluare

Rn(ξ ) = 2 ξ 2n+1

2n + 1+

ξ 2n+2

2n + 3+ . . .

< 2 · ξ 2n+1

2n + 1(1 + ξ 2 + ξ 4 + . . .) =

2 · ξ 2n+1

2n + 1· lim p→∞

ξ 2 p − 1

ξ 2 − 1=

ξ 2n+1

2n + 1· 2

1 − ξ 2

deci, folosind si (2.53)

Rn(ξ ) <ξ 2n+1

2n + 1

9

4, ξ ∈ (0, 1/3]. (2.54)

Astfel numarul ln z se va aproxima prin

ln z ≈ m · ln 2 − 2(u1 + . . . + un), (2.55)

unde

uk =ξ 2k−1

2k − 1, k = 1, . . . , n . (2.56)

Observatia 14 ˆ In (2.55) uk se calculeaz˘ a ca ın (2.56), iar pentru ln 2 se utilizeaz˘ a o valoare anterior calculat˘ a cu suficient de mare precizie (de ex.100 zecimale exacte).

Observatia 15 Din (2.54) si (2.56), deoarece ξ 2 ≤ 1/9 rezult˘ a

Rn(ξ ) ≤ 1

4un, (2.57)

relat ie care ne permite s˘ a evalu˘ am precizia de aproximare chiar ın timpul calculului expresiei din (2.55).

20



Capitolul 3

Aproximarea solutiilor

ecuatiilor si sistemelor de

ecuatii neliniare

Aproximarea solutiei (unice) ξ ∈ (a, b) a unei ecuatii neliniare deforma f (x) = 0, unde f : [a, b] → IR este o functie data, se realizeaza ınurmatoarele etape:

1. constructia unui sir (xn)n≥0 ⊂ [a, b] astfel ıncat xn → ξ ;

2. evaluarea erorii printr-o formula de tipul

|xn − ξ | ≤ E (n) → 0, n → ∞

Vom regasi aceste aspecte la metodele bisectiei, coardei, contractiilor sitangentei (a lui Newton) ce vor fi prezentate ın capitolul de fata.

3.1 Metoda bisectiei

Fie f : [a, b] → IR, continua astfel ıncat f (a)f (b) < 0. Presupunem,ın plus, ca radacina ξ a ecuatiei f (x) = 0 este unica ın (a, b). Consideram

urmatorul algoritm: fie I 0 = [a0, b0] = [a, b] si x0 =a0 + b0

2.

(a) Daca f (a)f (x) < 0, atunci b1 = x, a1 = a si luam I 1 = [a1, b1] = [a, x].(b) Daca f (x)f (b) < 0, atunci a1 = x, b1 = b si luam I 1 = [a1, b1] = [x, b].(c) Daca f (x) = 0 ınseamna ca ξ = x ; STOP.Repetand rationamentul, se obtin intervalele I 2 = [a2, b2], . . . , I n = [an, bn], . . .cu proprietatile

a0 ≤ a1 ≤ . . . ≤ an ≤ . . . (3.1)

21



b0

≥b1

≥. . .

≥bn

≥. . . (3.2)

f (an)f (bn) < 0, ∀ n ≥ 0, (3.3)

adica ξ ∈ I n = [an, bn] si

bn − an =b − a

2n. (3.4)

Teorema 5 (i) ˆ In condit iile anterioare, avem

∃ limn→∞ an = lim

n→∞ bn = ξ.

(ii) Dac˘ a definim ξ n =

an + bn

2 , ∀ n ≥ 0 atunci avem

|ξ n − ξ | ≤ b − a

2n, ∀ n ≥ 0.

Demonstratie. Din (3.1) si (3.2) si an ≤ bn rezulta

an ∈ (a0, bn), bn ∈ (an, b), ∀ n ≥ 0,

adica sirurile (an)n≥0 si (bn)n≥0 sunt marginite. Cum sirurile sunt si mono-tone rezulta ca ∃ lim

n

→∞

an = α si ∃ limn

→∞

bn = β . Cum an ≤ bn, ∀ n ≥ 0

obtinem α ≤ β . Vom arata ca α = β . Presupunem prin absurd ca α < β sifie 0 =

β − α

3. Din convergentele mentionate anterior rezulta pentru = 0

an → α ⇒ ∃ n1 ∈ IN astfel ıncat |an − α| < 0, ∀ n ≥ n1, (3.5)

bn → β ⇒ ∃ n2 ∈ IN astfel ıncat |bn − β | < 0, ∀ n ≥ n2, (3.6)

bn − an =b − a

2n→ 0 ⇒ ∃ n3 ∈ IN astfel ıncat bn − an < 0 (3.7)

pentru orice n ≥ n3. Atunci, pentru n ≥ n0 = max(n1, n2, n3) au loc

simultan relatiile (3.5) - (3.7). Atunci rezulta

−0 < α − an < 0, −0 < bn − β < 0, ∀ n ≥ n0.

Sumand inegalitatile din stanga, se obtine

−2

3(β − α) = −20 < bn − an − (β − α)

22



deciβ − α

3

< bn

−an,

∀n

≥n0

ceea ce este ın contradictie cu (3.7). Ramane atunci ca α = β . Apoi,trecand la limita ın (3.3) obtinem f (α)f (α) ≤ 0 adica f (α) = 0 si cum ξ era unica solutie din (a, b) obtinem α = ξ , deci concluzia de la punctul (i).Cea de-a doua parte a teoremei rezulta direct din consideratiile anterioare.

Observatia 16 Inegalitatea din teorema 5 (ii) permite o evaluare a erorii de aproximare. De exemplu, dac˘ a dorim s˘ a aproxim˘ am cu o eroare mai mic˘ a decˆ at 10−6 solut ia ecuat iei x = −ex folosind metoda bisect iei, proced˘ am astfel: consider am f (x) = x + ex; f (x) = 1 + ex > 0, ∀ x ∈ IR, deci f este strict cresc˘ atoare si cum f (0)f (−1) < 0 rezult˘ a ξ ∈ (−1, 0). Atunci

conform punctului (ii) avem c˘ a |ξ n − ξ | ≤ 12n ≤ 10−6 ⇒ n ≥ 20. Va fi deci suficient s˘ a se calculeze aproximat ia x20 prin metoda bisect iei pentru a avea eroarea dorit˘ a.

3.2 Metoda coardei

Fie f : [a, b] → IR de doua ori derivabila pe [a, b] cu

f (a)f (b) < 0 (3.8)

si

f (x) · f (x) = 0, ∀ x ∈ [a, b]. (3.9)

Lema 1 Avem urm˘ atoarele echivalent ¸e (i) f (a) · f (a) < 0 ⇔ f (b) · f (b) > 0;(ii) f (a) · f (a) > 0 ⇔ f (b) · f (b) < 0.

In raport cu proprietatea (3.9) si Lema 1, din cele 4 cazuri posibile lacapetele intervalului [a, b] raman doar doua. In raport cu aceste doua cazuriare loc urmatorul rezultat privind convergenta metodei coardei.

Teorema 6 ˆ In ipotezele (3.8)

−(3.9) si t inˆ and cont de Lema 1, unica solut ie

a ecuat iei f (x) = 0 poate fi obt inut˘ a ca limita sirului strict monoton din [a, b] definit astfel:(i) dac˘ a f (a)f (a) < 0 atunci

x0 = a

xn+1 = xn − f (xn)f (xn)−f (b) (xn − b), ∀ n ≥ 0

; (3.10)

23



Figura 3.1: Cazurile posibile pentru metoda coardei

(ii) dac˘ a f (b)f (b) < 0 atunci

x0 = b

xn+1 = xn − f (xn)f (xn)−f (a) (xn − a), ∀ n ≥ 0

. (3.11)

24



Teorema 7 Dac˘ a 0 < m1 ≤ |f (x)| , ∀ x ∈ [a, b], atunci avem evaluarea

|xn − ξ | ≤ |f (xn)|m1, ∀ n ≥ 1.

Demonstratie. Aplicand teorema lui Lagrange pe intervalul [xn, ξ ], obtinemun z ∈ (xn, ξ ) astfel ıncat f (ξ ) − f (xn) = (ξ − xn) · f (z). Atunci

|ξ − xn| =|f (xn)||f (z)| ≤ |f (xn)|

m1, ∀ n ≥ 1

si teorema este demonstrata.

Observatia 17 Formula din Teorema 7 permite o evaluare a posteriori a erorii de aproximare (dup˘ a calculul lui xn).

Observatia 18 Denumirea metodei coardei provine din interpretarea sa geometric˘ a (vezi de exemplu figura 3.2, pentru f (a) < 0, f (a) > 0). Avˆ and deja calculat˘ a iterat ia xn, adic˘ a cunoscˆ and punctul An(xn, f (xn)) urm˘ atoa-rea iterat ie xn+1 se obt ¸ine ca intersect ia dintre axa Ox si coarda A0An.Formula (3.10) (respectiv (3.11)) se obt ine scriind ecuat ia dreptei A0An si egalˆ and pe y cu 0.

Figura 3.2: Interpretarea geometrica a metodei coardei

Exemplul 9 S˘ a se aproximeze r˘ ad˘ acinile reale ale ecuat ¸iei x3 − 6x2 +10x = 4 utilizˆ and metoda coardei. Fie f (x) = x3 − 6x2 + 10x − 4. Cu

25



sirul lui Rolle, afl˘ am c˘ a r˘ ad˘ acinile ξ 1, ξ 2, ξ 3 ale lui f sunt toate reale si se g˘ asesc ın intervalele (0, 1), (1, 3) si respectiv (3,

∞). Vom exemplifica

pentru ξ 1 ∈ (0, 1). R˘ ad˘ acinile ecuat iei f (x) = 3x2 − 12x + 10 = 0 sunt x = 2 ±

23 ; atunci f (x) > 0, f (x) < 0, ∀ x ∈ [0, 1] si f (1)f (1) < 0.

Conform Teoremei 6 (ii) vom construi sirul x0 = 1,

xn+1 = xn − f (xn)

f (xn) − f (0)(xn − 0), ∀ n ≥ 0

cu xn → ξ 1. Cum f este strict descresc˘ atoare si pozitiv˘ a pe [0, 1], rezulta m1 = inf

x∈[0,1]|f ’ (x)| = f (1) = 1. Folosind formula anterioar˘ a calcul˘ am

primii trei termeni din sir prin

x1 = 1 − f (1)f (1) − f (0)

(1 − 0) = 0, 8

x2 = 0, 8 − f (0, 8)

f (0, 8) − f (0)(0, 8 − 8) ≈ 0, 68

x3 = 0, 68 − f (0, 68)

f (0, 68) − f (0)(0, 68 − 0) ≈ 0, 63

si cu evaluarea din teorema 7 rezult˘ a c˘ a

|ξ 1

−x3

| ≤ |f (0, 63)

|< 0, 15.

3.3 Metoda aproximatiilor succesive pe IR

Definitia 6 Fie g : [a, b] → IR. Un punct p ∈ [a, b] se numeste punct fixpentru g dac˘ a g( p) = p.Funct ia g se numeste α-contractie pe [a, b] dac˘ a ∃ α ∈ (0, 1) astfel ınc at

|g(x) − g(y)| ≤ α |x − y| , ∀ x, y ∈ [a, b]. (3.12)

Observatia 19 Dac˘ a g este o α-contract ie si are un punct fix, atunci el este unic. ˆ Intr-adev˘ ar, dac˘ a am avea p1, p2 ca puncte fixe pentru g ar rezulta

| p1 − p2| = |g( p1) − g( p2)| ≤ α | p1 − p2| ,

deci (cum α < 1) p1 = p2.

26



Observatia 20 Dac˘ a g : [a, b] → IR este derivabil˘ a si ∃ α ∈ [0, 1) astfel ınc at g(x) ≤ α, ∀ x ∈ [a, b] (3.13)

atunci g este o α-contract ie (se aplic˘ a teorema lui Lagrange pe un interval [x, y] ⊂ [a, b]).

Observatia 21 Dac˘ a g : [a, b] → [a, b] este continu˘ a, atunci g are un punct fix. ˆ Intr-adev˘ ar, fie h : [a, b] → IR, h(x) = g(x)−x. Cum h(a) = g(a)−a >0 si h(b) = g(b) − b < 0 rezult˘ a h(a) · h(b) ≤ 0 si cum h este continu˘ a,exist˘ a p ∈ [a, b] astfel ınc at h( p) = 0, adic˘ a g( p) = p.

Teorema 8 (Metoda aproximatiilor succesive)Fie g : [a, b] → [a, b] o α-contract ie, α ∈ [0, 1). Atunci ∀ x0 ∈ [a, b], sirul

(xn)n≥0 dat de xn = g(xn−1), ∀ n ≥ 1 (3.14)

converge c˘ atre unicul punct fix p al lui g. ˆ In plus, avem evaluarea

|xn − p| ≤ αn(b − a), ∀ n ≥ 1. (3.15)

Demonstratie. Avem succesiv

|xn − p| = |g(xn−1) − g( p)| ≤ α |xn−1 − p| ≤ . . . ≤

αn |x0 − p| ≤ αn(b − a), ∀ n ≥ 0.

Trecand la limita ın aceasta inegalitate, obtinem

limn→∞ |xn − p| ≤ (b − a) lim

n→∞αn = 0,

deci exista limn→∞xn = p si teorema este demonstrata.

Observatia 22 Rezultatul din Teorema 8 este valabil si pentru funct ii g : M → M , g α-contract ie pe M ⊆ IR, M mult ime ınchis a. Se poate ar˘ ata (vezi [ 26 ] si Teorema 11) c˘ a ın acest caz are loc evaluarea

|xn − p| ≤αn

1 − α |x0 − x1| , ∀ n ≥ 1. (3.16)

Exemplul 10 S˘ a se calculeze √

5 cu o eroare ≤ 10−2 utilizˆ and metoda aproximat iilor succesive.Fie g : (0, ∞) → (0, ∞), g(x) = 1

2 · (x + 5x ). Din tabelul de variat ¸ie al

funci ei g, rezult˘ a c˘ a √

5 este punct de minim global si, totodat˘ a punct fix

27



pentru g. Deci g((0, ∞)) ⊆ [√

5, ∞), dar g nu este o α-contract ie pe (0, ∞).Vom restrˆ ange ın acest caz domeniul lui g astfel ca g s˘ a fie α-contract ie pe

noul s˘ au domeniu. Se arat˘ a c˘ a g este o 12−contract ie pe [√ 5, ∞). Astfel definim g : [

√ 5, ∞) → [

√ 5, ∞) si din Observat ia 22 pentru x0 = 5/2

(de exemplu) si xn = g(xn−1) = . . . = gn(x0) obt inem c˘ a √

5 − xn

≤( 1

2 )n

1 − 12

x0 − 1

2(x0 +

5

x0)

=1

2n−1

x0

2− 5

2x0

=1

2n−1

1

4=

1

2n+1, ∀ n ≥ 1.

Cum 1

2n+1≤ 10−2 implic˘ a n ≥ 6, ınseamn a c˘ a

√ 5 este aproximat prin x6

cu o eroare mai mic˘ a decˆ at 10−2, unde x6 = g(x5) = . . . = g6(x0).

3.4 Metoda lui Newton pe IR

Fie f : [a, b] → IR de clasa C 2([a, b]). Presupunem ca ecuatia f (x) = 0are o unica solutie ξ ∈ [a, b].

Figura 3.3: Interpretarea geometrica a metodei lui Newton

Din interpretarea geometrica a metodei (xn este intersectia tangenteila graficul lui f ın (xn−1, f (xn−1)) cu axa Ox - vezi figura 3.3) se deduceformula de recurenta

xn

= xn−1 −

f (xn−1)

f (xn−1),

∀n

≥1. (3.17)

Teorema 9 (Metoda Newton)(i) ˆ In condit iile anterioare, dac˘ a ın plus

f (ξ ) = 0 (3.18)

28



atunci ∃ δ > 0 astfel ınc at pentru orice x0 ∈ [ξ − δ, ξ + δ ], sirul (xn)n≥1 dat de (3.17), converge la ξ .

(ii) Dac˘ a m1, M 2 > 0 sunt date de

m1 = inf x∈[a,b]

f (x) > 0, M 2 = sup

x∈[a,b]

f (x) (3.19)

atunci avem evalu˘ arile

|ξ − xn| ≤ M 22m1

(xn − xn−1)2, ∀ n ≥ 1, (3.20)

|ξ − xn| ≤ M 22m1

(ξ − xn−1)2, ∀ n ≥ 1. (3.21)

Observatia 23 Un dezavantaj al acestei metode este alegerea aproximat iei init iale x0 care trebuie s˘ a satisfac˘ a condit ia xo ∈ [ξ − δ, ξ + δ ] , valorile ξ si δ nefiind cunoscute a priori. Acest neajuns este eliminat de urm˘ atorul rezultat.

Teorema 10 (Metoda tangentei)Fie f : [a, b] → IR, de dou˘ a ori derivabil˘ a astfel ınc at f (a)f (b) < 0, f (x) =0 si f (x) = 0, ∀ x ∈ [a, b]. Atunci sirul (xn)n≥0 dat de (3.17) cu x0 ∈ [a, b]astfel ınc at

f (x0)f (x0) > 0,

converge (strict monoton) la ξ , unica solut ie a ecuat iei f (x) = 0.

Exemplul 11 S˘ a se aproximeze solut iile reale ale ecuat iei ex + 2x + 1 = 0utilizˆ and metoda tangentei. Fie f : IR → IR, f (x) = ex + 2x + 1. Cu sirul lui Rolle afl˘ am c˘ a ecuat ia are o singur˘ a r˘ ad˘ acin˘ a real˘ a x∗ ∈ (−1, 0). Cum f (x) > 0, f (x) > 0, ∀ x ∈ [0, 1], se construieste sirul

x0 = 0, xn = xn−1 − f (xn−1)

f (xn−1)→ x∗.

Se observ˘ a c˘ a f (0)f (0) > 0. Obt inem

x1 = 0 −f (0)

f (0) = −2

3 ≈ −0.6666666,

x2 = −2

3− f ( 2

3 )

f ( 23 )

≈ −0, 7383188,

m1 = inf x∈[−1,0]

f (x) =

f (−1) = e−1 + 2 ≈ 2, 3678905

29



M 2 = supx∈[−1,0]

f (x)

=f (0)

si conform (3.20) avem

|x∗ − x2| ≤ m2

2m2|x2 − x1|2 ≈ 0, 001084.

3.5 Metoda aproximatiilor succesive si metoda

Newton pentru sisteme de ecuatii neliniare

In aceasta sectiune vom prezenta extinderi ale metodelor aproximatiilorsuccesive si Newton pentru aproximarea solutiilor sistemelor de ecuatiineliniare de forma F : D ⊂ IRn → IRn, cu n ≥ 2.

Definitia 7 Fie α ∈ (0, 1) fixat, D ⊂ IRn domeniu, o aplicat ie

g : D ⊆ IRn

→ IRn

se numeste α-contractie dac˘ a g(y) − g(x) ≤ α · x − y , ∀ x, y ∈ D

(unde · este o norm˘ a pe IRn). Un element p ∈ D se numeste punct fixal lui g dac˘ a g( p) = p.

Teorema 11 (Metoda aproximatiilor succesive pe IRn)Fie D ⊆ Rn, D domeniu ınchis, α ∈ [0, 1), g : D → IRn o α-contract ie si x0 ∈ D. Definim xk+1 = g(xk), ∀ k ≥ 0 si presupunem c˘ a xk ∈ D, ∀ k ≥ 0.Atunci avem (i) Sirul (xk)k≥0 converge la unicul punct fix p al lui g ın D;(ii) Avem evaluarea

xk − p ≤ αk

1 − αx1 − x0 , ∀ n ≥ 1. (3.22)

Teorema 12 Fie D ⊆ IRn compact˘ a, g : D → IRn o α-contract ie astfel ınc at x0 ∈ D, x1 = g(x0) ∈ D si

dist((x1, F r D) ≥ dα

1 − α,

unde d este diametrul lui D. Atunci, pentru orice k ≥ 0, xk ∈ D.

Observatia 24 Dac˘ a D

⊆IRn ınchis astfel ınc at

∃V o vecin˘ atate a lui D,

g : D → IRn cu g ∈ C 1(V ) si ∃ α ∈ (0, 1) astfel ınc at g(x) ≤ α ∀ x ∈ D,atunci g este o α-contract ie. ˆ Intr-adev˘ ar, aplicˆ and teorema cresterilor finite (vezi [26 ]) avem c˘ a

g(x) − g(y) ≤ x − y supξ∈[x,y]

g(ξ ) ≤ α · x − y , ∀ x, y ∈ V.

30



Exemplul 12 . Fie funct ¸ia g : [1, 2] × [0, 1] → IR2, dat˘ a de g(x, y) =

(1 +3x2 + y2

20,

2x2 + 5y2

20). S˘ a se arate c˘ a g este α-contract ie ın raport cu

·∞. Obt ¸inem c˘ a

g(x, y) =

3x

10

y

10

x

5

y

2

.

Cum toate cele patru funct ii care apar sunt pozitiv cresc˘ atoare obt inem

g(x, y)∞ = max

(x,y)∈D3x

10

+ y

10

, x

15

+y

2

= max3 · 2

10+

1

10,

2

5+

1

2 =

max7

10 ,

9

10 =

9

10 ∈ [0, 1).

si conform Observat iei 24, g este o9

10-contract ie pe [1, 2] × [0, 1].

Vom da fara demonstratie urmatorul rezultat (pentru detalii vezi [6] si [22]).

Teorema 13 (Metoda Newton-Raphson)Fie F = (F 1, . . . , F n) : D ⊆ IRn → IRn, x0 ∈ D , h > 0 astfel ınc at B[x0, h] ⊂ D. Dac˘ a sunt ındeplinite condit iile:

(i) ∃ F (x0)−1

= Γ0 si Γ0 ≤ A0;(ii) ∃ B0 > 0 astfel ınc at

Γ0F (x0)

≤ B0 ≤ h

2 ;

(iii) ∃ C 0 > 0 astfel ınc at

max1≤k≤n

max1≤i≤n

n j=1

∂ 2F k(x)

∂xi∂x j

≤ C 0, ∀ x ∈ B[x0, h]

(iv) µ0 = 2nA0B0C 0 < 1.Atunci sirul definit de metoda Newton xk+1 = xk − (F (xk))−1 · F (xk),k ≥ 0 este bine definit si converge la un element ξ ∈ B[x0, h] care este solut ie a sistemului de ecuat ii F (x) = 0.

Observatia 25 Spre deosebire de metoda aproximat iilor succesive, metoda Newton nu necesit˘ a nici o transformare a sistemului init ial. ˆ In plus, viteza

de convergent a a metodei Newton este mult mai bun˘ a decˆ at cea dat˘ a de metoda aproximat iilor succesive (convergent a p˘ atratic˘ a fat a de convergent a liniar˘ a; vezi pentru detalii [ 12 ]).

Observatia 26 L˘ asˆ and la o parte determinarea aproximat iei x0 si a bilei B[x0, h] pe care sunt date condit iile ce asigur˘ a convergent a, metoda Newton

31



prezint˘ a unele dezavantaje. Astfel, la fiecare itert ie trebuie calculat i F (xk)si F

(xk) si, de asemenea, trebuie rezolvat sistemul F

(xk)

·y = F (xk), ceea

ce ınseamn˘ a efort computat ional mare si poate duce la acumularea de erori.Aceste inconveniente pot fi ınl˘ aturate folosind metoda Newton simplificat˘ a (modificat˘ a) pe care o prezent˘ am ın continuare.

Teorema 14 ˆ In condit iile Teoremei 13 , sirul definit de

xk+1 = xk − (F (x0))−1F (xk), k ≥ 0

converge la un element ξ ∈ B[x0, h] care este solut ie a sistemului de ecuat ii F (x) = 0.

Exemplul 13 . Fie sistemul x − x2

20− y

20= 2

y − y2

20− x

20− 1 = 0

si

D = (x, y) ∈ IR2, (x − 5

2)2 + (y − 3

2)2 < 2.

S˘ a se verifice dac˘ a sunt ındeplinite condit iile din teorema de convergent a de la metoda Newton relativ la aproximat ia init ial˘ a (x0, y0) = ( 5

2 , 32 ).

Fie F : D

→IR2, F (x, y) = (x

−x2

20 −y

20 −2, y

−y2

20 −x

20 −1). Avem

succesiv (ın raport cu · ∞ vectorial˘ a si matriceal˘ a):

(i) calcul˘ am Γ0 = (F (x0, y0))−1 si afl am A0 =180

127;

(ii) calcul˘ am apoi Γ0F (x0, y0) si afl am B0 =81

254;

(iii) calcul˘ am C si afl am C =1

10.

ˆ In plus, pentru constanta µ0 avem

µ0 = 2nA0B0C 0 <1

5< 1 ⇒

deci, se poate aplica Teorema 13.

32



3.6 Exercitii

1. Aproximati radacinile ecuatiei √ x = cos x pe [0, 1] prin x3 dat demetoda bisectiei.

2. Aproximati radacinile ecuatiei x3 − 7x2 + 14x − 6 = 0 prin metodabisectiei cu o eroare mai mica decat 10−2 pe intervalele(a) [0, 1], (b) [1, 3.2], respectiv (c) [3.2, 4].

3. Aproximati 3√

25 prin metoda bisectiei cu o eroare mai mica decat10−4.

4. Sa se discute dupa parametrul real m numarul de radacini reale pen-tru ecuatia

|ln x

| −m

√ x = 0. Pentru m = 2, aproximati radacina

subunitara folosind metoda bisectiei cu o eroare mai mica de 10−3.

5. Aproximati 3√

25 prin metoda aproximatiilor succesive cu o eroare maimica decat 10−4. Comparati rezultatul cu cel obtinut la exercitiul 3.

6. Aproximati radacinile urmatoarelor ecuatii prin metoda aproximatiilorsuccesive cu o eroare mai mica de 10−5

(a) 3x2 − ex = 0;

(b) x − cos x = 0;

(c) x =

5

x2 + 2;

(d) x =2 − ex + x2

3;

(e) x2 + 10 cos x = 0.

7. Functia f (x) = tgπx − 6 are ca radacina pe1

πarctg6 0.447431543.

Fie p0 = 0 si p1 = 0.48. Aproximati aceste radacini prin x10 dat de

(a) metoda coardei;

(b) metoda bisectiei.

Prin care din cele doua metode se obtine o eroare mai mica?

8. Sistemul neliniar x2 − 10x + y2 + 8 = 0xy2 + x − 10y + 8 = 0

33



poate fi adus ın forma

x = g1(x, y) = x2

+ y2

+ 810

y = g2(x, y) =xy2 + x + 8

10

unde g = (g1, g2) : D → IR, cu D = (x, y) ∈ IR2, 0 ≤ x, y ≤ 1.5.

(a) Aratati ca g are un unic punct fix pe D.Indicat ie. Se va arata mai ıntai ca g este o α- contractie pe Dın raport cu · ∞, apoi ca g(D) ⊆ D.

(b) Determinati numarul minim de iteratii necesare astfel ca eroareaaproximarii prin metoda contractiilor sa fie mai mica decat 10−6.

34



Capitolul 4

Aproximare, interpolare si

derivare numerica

Se vor analiza ın cadrul acestui capitol urmatoarele doua probleme impor-tante ale analizei numerice.

Problema 1 Dat˘ a fiind o funct ie f : [a, b] → IR cu anumite propriet˘ at i,s˘ a se determine un polinom P ∈ IR[X ] astfel ınc at

f − P ≤ ε, (4.1)

unde · este o norm˘ a pe spat iul vectorial al funct iilor f : [a, b] → IR, iar

ε e o precizie dat˘ a.

Problema 2 Dat˘ a fiind o funct ie f : [a, b] → IR, diviziunea : (a = x1 <x2 < ... < xn = b), s˘ a se determine o funct ie ϕ : [a, b] → IR polinomial˘ a sau polinomial˘ a pe buc˘ at i astfel ınc at

f (xi) = ϕ(xi), i = 1,...,n (4.2)

si f − ϕ ≤ ε. (4.3)

Problema 1 se refera la aproximarea unei functii prin functii polinomiale,

iar cea de-a doua la interpolarea si aproximarea unor functii prin polinoamesau functii ce sunt polinomiale pe intervale. Se folosesc functii polinomi-ale sau polinomiale pe bucati pentru ca valoarea functiei ıntr-un punct sepoate calcula foarte simplu (vezi Capitolul 2, Schema lui Horner). In cazulProblemei 2, functia f este de cele mai multe ori cunoscuta doar ın punctelediviziunii .

35



4.1 Aproximarea cu polinoame Bernstein

Incepem aceasta sectiune cu un rezultat auxiliar.

Lema 2 Pentru ∀ x ∈ IR si n ∈ IN sunt adev˘ arate urm˘ atoarele egalit˘ at i.

(i)n

k=0

kC knxk (1 − x)n−k = nx;

(ii)n

k=0

k2C knxk (1 − x)n−k = nx + n (n − 1) x2;

(iii)n

k=0(nx − k)2 C knxk (1 − x)n−k = nx (1 − x).

Teorema 15 Fie f : [0, 1] → IR, continu˘ a si (Bn)n≥1 un sir de funct ii polinomiale definite prin

Bn (x) =n

k+0

f k

n

C knxk (1 − x)n−k , ∀ x ∈ [0, 1]. (4.4)

Atunci sirul (Bn)n≥1 converge uniform la f pe [0, 1].

Observatia 27 Polinomul Bn(x) se numeste polinom Bernstein de or-din n asociat funct iei f pe intervalul [0, 1].

Extinderea aproximarii date de Teorema 15 la un interval arbitrar [a, b] esteprezentata ın rezultatul urmator.

Teorema 16 Fie f : [a, b]→

IR continu˘ a. Atunci exist˘ a (P n

)n≥1

un sir de polinoame cu coeficient i reali astfel ıncˆ at f − P n∞ −→ 0, n → ∞.

Demonstratie. Fie (Bn)n≥1 polinoamele Bernstein din Teorema 15 aso-ciate functiei F : [0, 1] → IR, F (t) = f (a + t (b − a)) si

P n (x) = Bn

x − a

b − a

, x ∈ [a, b].

Folosind Teorema 15 obtinem succesiv

f − Bn∞ = supx∈[a,b]

| f (x) − P n (x) |= supt∈[0,1]

| F (t) − Bn (t) |−→ 0, n → ∞

ceea ce ıncheie demonstratia teoremei.

Observatia 28 Teorema 16 poate fi extins˘ a ın modul urmator (pentru demonstrat ¸ii vezi [ 6 ], [ 22 ]): dac˘ a f : [a, b] → IR, f ∈ C k ([a, b]), atunci

∃ (P n)n≥1 un sir de polinoame astfel ıncˆ at f (i) − P (i)n ∞ → 0, pentru

∀ i = 0,...,k.

36



Observatia 29 Asa cum rezult a din demonstrat ia Teoremei 15, convergen-t a sirului (Bn)n≥1 este foarte lent˘ a, ceea ce reprezint˘ a un dezavantaj pen-

tru aproximarea cu polinoame Bernstein. ˆ In plus, nu avem posibilitatea evalu˘ arii a priori a diferent ei | f − Bn |. Dar, avantajul metodei const˘ a ın faptul c˘ a ∀ x ∈ [a, b], Bn(x) este o combinat ¸ie convex˘ a a valorilor f

kn

, k = 0,...,n, proprietate ce p˘ astreaz˘ a ”alura” (shape) funct iei f (adic˘ a

intervalele de convexitate si concavitate).

Observatia 30 Cum Bn

kn

= f

kn

, rezult˘ a c˘ a polinomul Bernstein nu

are proprietatea de interpolare (4.2), ci doar pe cea de aproximare dat˘ a de Teorema 15.

4.2 Interpolare cu polinoame Lagrange

In aceasta sectiune ne vom referi la Problema 2. Astfel, fiind date valorile luif ın n puncte distincte xi, i = 1,...,n, se pune problema aproximarii printr-un polinom al carui grafic sa ”treaca” prin punctele Ai (xi, f (xi)) , i =1,...,n (vezi figura 4.1). Urmatorul rezultat ne asigura ca exista un astfelde polinom.

Figura 4.1: Interpolare polinomiala

Teorema 17 Fie f : A ⊂ IR → IR si x1,...,xn ∈ A puncte distincte.Atunci exist˘ a si este unic determinat polinomul P ∈ IR[X ] de grad cel mult n − 1 astfel ınc at

f (xi) = P (xi) , ∀ i = 1,...,n. (4.5)

37



Polinomul P poate fi reprezentat sub forma

P (x) =

n j=1

f (xi)

ni=1

x − xix j − xi

, x ∈ IR. (4.6)

Definitia 8 Polinomul

Ln (f ; x) =n

j=1

f (x j) Ln,j (x) . (4.7)

se numeste polinomul de interpolare al lui Lagrange asociat funct iei f si punctelor xi.

Fie ω : [a, b] → IR polinomul unitar de grad n definit de dat de

ω (x) = (x − x1) ... (x − xn) . (4.8)

In punctele xi, i = 1, . . . , n , Ln(f ; x) are proprietatea de interpolare (4.2);ın celelalte puncte x ∈ [a, b] are loc urmatorul rezultat (ce asigura ın anu-mite conditii si proprietatea de aproximare).

Teorema 18 Dac˘ a f ∈ C n ([a, b]), atunci ∀ x ∈ [a, b], x = xi, pentru ∀ i = 1,...,n, exist˘ a ξ x ∈ [a, b] astfel ınc at

f (x) − Ln (f ; x) =

f (n) (ξ x)

n! ω (x) . (4.9)

Lema 3 Dac˘ a h = maxi=2,...,n

| xi − xi−1 |> 0, atunci

ω∞ ≤ hn

4(n − 1)!

Demonstratie. Fie x ∈ [a, b] fixat arbitrar si i ∈ 1, . . . , n astfel ıncatx

∈(xi

−1, xi) . Atunci avem

| (x − xi−1) (x − xi) |≤ (xi − xi−1)2

4≤ h2

4. (4.10)

Cum ınsa, pentru j ∈ 1, . . . , i − 2 avem

| x − x j |= x − x j ≤ (i − j) h, (4.11)

38



iar pentru j ∈ i + 1, . . . , n

| x − x j |= x − x j ≤ ( j − i + 1) h, (4.12)

din (4.10)-(4.12) obtinem

| ω (x) |≤ hi−2 (i − 1)!h2

4hn−i (n − i)! ≤ (n − 1)!

hn

4

si teorema este demonstrata.Din Teorema 18 si Lema 3 obtinem urmatorul rezultat.

Teorema 19 ˆ In condit iile din Teorema 18, are loc relat ia

f

−Ln (f ;

·)∞ ≤

f (n)∞4n

hn. (4.13)

Observatia 31 Dac˘ a f ∈ C ∞ ([a, b]), f (n)∞ ≤ M, ∀ n ≥ 1 si h < 1atunci din (4.13) rezult˘ a c˘ a

f − Ln (f ; ·) ∞ → 0, n → ∞

adic˘ a sirul polinoamelor Lagrange converge uniform la f , ceea ce ınseamna ca aceste polinoame au si proprietatea de aproximare (4.1).

Observatia 32 Un mare dezavantaj al aproxim˘ arii cu polinoame Bern-stein sau Lagrange este faptul c˘ a gradul polinomului creste odat˘ a cu num˘ arul de puncte din [a, b]. Acest lucru face ca valoarea polinomului ıntr-un punct s˘ a se calculeze cu mult efort si erori de calcul mari. Acest inconvenient poate fi eliminat prin utilizarea funct iilor polinomiale pe bucat i. Un exem-plu ın acest sens este prezentat ın sect iunea urm˘ atoare.

4.3 Functii spline cubice

Definitia 9 Se numeste functie spline cubica o funct ie g de clas˘ a C 2 ([a, b]), polinomial˘ a de grad 3 pe buc˘ at i, cu urm˘ atoarele propiet˘ at i

g|[xi−1,xi](x) = gi (x) =

3l=0

ail (xi − x)

l

, i = 1,...,n, (4.14)

g (xi) = f (xi) ; (4.15)

g

(xi) = f

(xi) ; (4.16)

39



g

(xi) = f

(xi) ; (4.17)

g (a) = g (b) = 0, (4.18)

unde a = x0 < x1 < · · · < xn = b este o diviziune a intervalului [a, b].

Observatia 33 O funct ie g ce verific˘ a relat iile (4.14) − (4.18) se numeste functie spline cubica libera, iar dac˘ a ın loc de relat ia (4.18) se impun condit iile

g (a) = f (a) , g (b) = f (b) , (4.19)

g se numeste functie spline cubica fixata.

Urmatorul rezultat asigura, ın conditiile date existenta unei functii splinecubice (libere).

Teorema 20 Exist˘ a si este unic˘ a o funct ie spline cubic˘ a liber˘ a.

Observatia 34 Un rezultat analog se demonstreaz˘ a pentru existent a unei funct ii spline cubice fixate (vezi (4.19)). Pentru detalii si extinderi ale not iunii de funct ie spline se poate consulta lucrarea [ 23 ].

4.4 Derivarea numerica

Necesitatea derivarii numerice apare direct legata de situatiile ın care functia

f : [a, b] → IR are o expresie prea complicata sau cand valorile ei suntcunoscute doar ın punctele x0, x1,...,xn ale unei diviziuni ∆ a intervalului[a, b]. Dar operatorii de derivare numerica apar si ın constructiile schemelorcu diferente finite utilizate ın aproximarea solutiilor ecuatiilor cu derivatepartiale (vezi [11]). Prin derivare numerica ıntelegem determinarea uneiformule de calcul care sa aproximeze valoarea derivatei lui f ıntr-un punctdin intervalul [a, b] si evaluarea erorii de aproximare. Asemenea formule seobtin de obicei plecand de la o functie de interpolare asociata lui f si divi-ziunii ∆. Vom prezenta ın aceasta sectiune o varianta care utilizeaza inter-polarea Lagrange. Fie deci f : [a, b] → IR si ∆ : a = x0 < x1 < ... < xn = bo diviziune echidistanta a lui [a, b],

xi+1 − xi = h =b − a

ni = 1,...,n − 1, (4.20)

x ∈ [xi, xi+1] ⊂ [a, b] arbitrar fixat si q ≥ 0 dat de

q =x − x0

h. (4.21)

40



Daca ω : [a, b] → IR este polinomul

ω (x) = (x − x0) (x − x1) ... (x − xn) (4.22)atunci avem

ω (x) = hn+1q (q − 1) ... (q − n) , (4.23)

ω

(xi) = (−1)n−i hn (i)! (n − i)! (4.24)

Sa notamα (q, n) = q (q − 1) (q − 2) ... (q − n) (4.25)

si fie Ln+1 (x) polinomul Lagrange de interpolare

Ln+1

(x) =n

i=0

f (xi)

(x − x0) ... (x − xi−1) (x − xi+1) ... (x − xn)

(xi − x0) ... (xi − xi−1) (xi − xi+1) ... (xi − xn).

(4.26)Utilizand notatiile si formulele anterioare, un calcul simplu ne arata ca

Ln+1 (x) =

ni=0

(−1)n−1 f (xi)

i! (n − i)!· α (q, n)

q − i. (4.27)

In plus, din (4.21) rezulta

x = x (q ) = x0 + hq, xq =dx

dq = h. (4.28)

Fie Ln+1 (q ) polinomul ın variabila q dat de

Ln+1 (q ) = Ln+1 (x (q )) . (4.29)

Din (4.28) si (4.29) obtinem

d

dx(Ln+1 (x)) =

d

dq

1

hLn+1 (q )

=

1

h

ni=0

(−1)n−i f (xi)

i! (n − i)!

d

dq

α (q, n)

q − i

. (4.30)

Fie Rn+1 (x) dat de

Rn+1 (x) = f (x) − Ln+1 (x) . (4.31)

Din teorema 18 rezulta ca daca f este de clasa C n+1([a, b]) exista ξ x ∈ (a, b)astfel ıncat

Rn+1 (x) =f n+1 (ξ (x))

(n + 1)!ω (x) . (4.32)

41



Daca f este de clasa C n+2 pe [a, b] derivand ın raport cu x relatia siınlocuind x cu xi rezulta

R

n+1 (xi) =1

(n + 1)!ω

(xi) f n+1 (ξ ) =

(−1)n−i hni! (n − i)!

(n + 1)!f n+1 (ξ x) (4.33)

Din (4.31) obtinem atunci

f

(xi) = L

n+1 (xi) + R

n+1 (xi) , i = 0,...,n. (4.34)

Partea dreapta a egalitatii 4.34 reprezinta o formula de derivare numerica

pentru f ın punctele diviziunii ∆ ın care R

n+1 (xi) reprezinta restul deaproximare. De exemplu, pentru n = 2 obtinem

L3 (x) =1

2f (x0) (q − 1) (q − 2) − f (x1) q (q − 2) +

1

2f (x2) q (q − 1)

L

(x) =1

h[1

2f (x0) (2q − 3) − f (x1) (2q − q )]

si resturile de aproximare de ordin 2

R

3 (x0) =1

3h2f (3) (ξ 0) , R

3 (x1) = −1

6h2f (3) (ξ 0) , R

3 (x2) =1

3h2f (3) (ξ 2) .

Decif

(x0) 1

2h((−3) f (x0) + 4f (x1) − f (x2)) ,

f

(x1) 1

2h(−f (x0) + f (x2)) ,

f

(x2) 1

2h(f (x0) − 4f (x1) + 3f (x2)) ,

toate aproximarile fiind de ordin 2.

Observatia 35 Modul ın care aplic˘ am rezultatele anterioare este urm˘ ato-rul. Fiind dat x ∈ (a, b) arbitrar consider˘ am un interval [α, β ] ⊂ (a, b) ce cont ine pe x. Consider˘ am apoi o diviziune a acestui interval α = x

0< x

1<

... < xN = β astfel ınc at x s˘ a fie unul din punctele diviziunii. Aplic˘ am apoi rat ¸ionamentele anterioare pentru [α, β ] si diviziunea considerat˘ a. De obicei lu˘ am

N = 2, α = x0 = x − ε, β = x2 = x + ε, x1 = x,

unde ε > 0 este ales ın raport cu precizia dorit˘ a.

42



4.5 Exercitii

1. Folosind polinomul de interpolare Lagrange aproximati:

(a) f (8.2) daca se dau nodurile de interpolare x0 = 8.1,x1 = 8.3, x2 = 8.4, x3 = 8.7 si f (x) = (x + 1)ex.

(b) f (0.8) pentru x0 = 0.6, x1 = 0.7 x2 = 0.9, x3 = 1.0 sif (x) = sin(x2 − 2).

2. Comparati rezultatul obtinut la exercitiul 1 cu valoarea exacta.

3. Comparati diferentele ıntre rezultatul obtinut la exercitiul 1 si rezul-tatul obtinut la exercitiul 2 cu eroarea data de teorema 17.

4. Aproximati f (1.12) prin polinomul Lagrange corespunzator dacaf (1.0 ) = 0.1924, f (1.05) = 0.2414, f (1.10) = 0.2933, f (1.15) =0.3492.

5. Sa se construiasca polinomul Lagrange ce interpoleaza functia f (x) =x5 ın nodurile x0, x1, si x4, unde xi = i

6 , i = 0, . . . , 6.

Indicat ie. P 0,1,4(x) = (x−x1)(x−x4)(x0−x1)(x0−x4) f (x0) + (x−x0)(x−x4)

(x1−x0)(x1−x4) f (x1) +(x−x0)(x−x1)

(x4−x0)(x4−x1) f (x4)

Observatia 36 O alt˘ a metod˘ a de a construi polinomul de interpolare (Lagrange) este metoda diferentelor finite a lui Newton. Astfel,P

n(x) poate fi scris ca:

P n(x) = f [x0] +

nk=1

f [x0, x1,...,xk](x − x0)...(x − xk−1)

unde f [xi] = f (xi) este diferent a finit˘ a de ordin 0,

f [xi, x j ] =f [x j] − f [xi]

x j − xi, i < j este diferent a finit˘ a de ordinul 1

si, ın general

f [xi, xi+1,...,xi+ j ] =f [xi+1, xi+2,...,xi+ j] − f [xi, xi+1,...,xi+ j−1]

xi+ j − xi

este diferent a finit˘ a de ordinul j.

6. Folosind polinoamele de interpolare Newton descrise mai sus, aproximatif −1

3

daca se dau f (−1.0) = −1, f (−0.75) = −0.421875, f (−0.5) =

−0.125 si f (0.25) = 0.015625

43



7. Daca P n(x) = f [x0] + f [x0, x1](x − x0) + a2(x − x0)(x − x1)++a3(x

−x0)(x

−x1)(x

−x2) + ... + an(x

−x0)(x

−x1)...(x

−xm−1)

aratati ca a2 = f [x0, x1, x2]. Cat este a3?

8. Determinati functia spline cubica fixata cu proprietatile f (0) = 0,f (1) = 1, f (2) = 2 si S (0) = S (2) = 1.

9. Daca punctele (xi, f (xi))i=1,n sunt asezate pe o dreapta, ce sepoate spune despre functia spline cubica (libera sau fixata) atasatafunctiei f ?Indicat ie. Utilizati exercitiile ?? si 8.

10. Fie f : [a, b] → IR si nodurile x0 = a < x1 < x2 = b. O functie splinepatratica este o functie S cu proprietatea:

S (x) = a0 + b0(x − x0) + c0(x − x0)2, x ∈ [x0, x1]a1 + b1(x − x1) + c1(x − x1)2, x ∈ [x1, x2]

astfel ıncat S (xi) = f (xi), i = 0, 2, S ∈ C 1[x0, x2] Sunt suficienteaceste conditii pentru a determina functia S ? Ce se poate spunedespre solutia obtinuta daca cerem ca S ∈ C 2[x0, x2]?

44



Capitolul 5

Integrare numerica

Fie f : [a, b] → IR o functie continua. Pentru a calcula integrala definita ba

f (x)dx se pot utiliza mai multe metode, ca de exemplu

1. se foloseste definitia integralei cu suma Riemann sau cu sumele Dar-boux;

2. daca se cunoaste o primitiva F a lui f , se foloseste formula Leibnitz-

Newton

ba

f (x)dx = F (b) − F (a);

3. se dezvolta f ın serie Taylor ın jurul punctului x0 ∈ [a, b], adica

f (x) =∞

n=0

an(x − x0)n,

iar apoi se integreaza seria data termen cu termen pe [a, b].

Dar, de multe ori, aceste metode sunt sau foarte laborioase sau chiar im-posibil de aplicat, deci impracticabile chiar si pentru unele functii simplecum ar fi

f : [0, 1] → IR; f (x) = e−x2; g : [0, π] → IR; g(x) = sin(x2).

In aceste cazuri, se aproximeaza valoarea integralelor cu o suma de tipul

ni=1

αif (xi) (5.1)

care se numeste formula de cuadratura.

45



5.1 Formule de cuadratura de tip Newton-Cotes

5.1.1 Aproximarea integralelor pe IR

Daca f ∈ C n([a, b]), fie t1, t2,...,tn ∈ [−1, 1] fixati si Ln polinomul deinterpolare al lui Lagrange asociat lui f si punctelor

x j =b + a

2+

b − a

2t j, j = 1,..,n (5.2)

Cum Ln aproximeaza pe f ,

ba

Ln(x)dx va aproxima

ba

f (x)dx. Eroarea

acestei aproximari este data de

rn(f ) = ba

(f (x) − Ln(x))dx.

Din Teorema 18 avem ca

|f (x) − Ln(x)| ≤ f (n) ∞n!

|ωn(x)|, ∀ x ∈ [a, b],

de unde, folosind formula anterioara rezulta ca

|rn(f )| ≤ f (n) ∞n!

ba

|ωn(x)|dx.

Cu schimbarea de variabila

x =b + a

2+

b − a

2t, t ∈ [−1, 1] (5.3)

rezulta apoi ba

|ωn(x)|dx =

1

−1

ωn

b + a

2+

b − a

2t

b − a

2dt =

1

−1

n

j=1 b + a

2+

b − a

2t − x j

b − a

2dt =

b − a

2 n+1

1

−1

n

j=1

|t − tj|dt

(s-a folosit relatia (5.2)). Asadar

|rn(f )| ≤ f (n) ∞n!

b − a

2

n+11

−1

n j=1

|t − t j|dt. (5.4)

46



Presupunand ca toate puctele t j, j = 1,..,n sunt distincte rezulta ca sipunctele x j, j = 1,..,n sunt distincte, iar din (4.7) obtinem

Ln(x) =n

j=1

f (x j)n

i=1,i= j

x − xi

x j − xi

si

b a

Ln(x)dx =

n j=1

1 −1

f

b + a

2+

b − a

2t j

b − a

2

ni=1,i= j

t − ti

t j − tidt =

=b−

a

2

n j=1

f b + a

2 +b−

a

2 t j1

−1

ni=1,i= j

t−

ti

t j − tidt. (5.5)

Decib

af (x)dx este aproximata de integrala (5.5), cu o eroare data de (5.4).

Teorema 21 (Formula dreptunghiurilor) Fie f:[a,b] → IR, f ∈ C ([a, b])si ∆ : (a = x0 < ... < xn = b) o diviziune echidistant˘ a a intervalului [a, b].Definim num˘ arul σ1(f, n) prin

σ1(f, n) =b − a

n

n

i=0

f xi + xi+1

2 . (5.6)

Atunci σ1(f, n) aproximeaz˘ a b

a f (x)dx cu o eroare majorat˘ a de

(b − a)2

4n f ∞ . (5.7)

Demonstratie. In relatia (5.5), consideram n = 1 si t1 = 0. Obtinem ca

b

a

f (x)dx b − a

2f

b + a

2 1

−1

1dt = (b − a)f

b + a

2 ,

cu o eroare majorata de

f ∞1!

b − a

2

21

−1

|t|dt =(b − a)2

4f ∞ (5.8)

47



Asadar aplicand (5.8) pentru intervalul [xi, xi+1] vom avea

xi+1 xi

f (x)dx (xi+1 − xi)f xi + xi+1

2 = b − a

n· f xi + xi+1

2 ,


(xi+1 − xi)2

4sup

x∈[xi,xi+1]| f (x) |≤ (b − a)2

4n2f ∞,

unde

xi = a +b − a

ni, ∀ i = 0,...,n.

Deci b a

f (x)dx =

n−1i=0

xi+1 xi

f (x)dx b − a

n

n−1i=0

f

xi + xi+1

2

,


n · (b − a)2

4n2f ∞ =

(b − a)2

4nf ∞

si teorema este complet demonstrata.

Observatia 37 Formula (5.6) este exact˘ a pentru f ∞ = 0 adic˘ a pentru funct iile constante.

Figura 5.1: Metoda dreptunghiurilor

48



Observatia 38 Relat ia

b a

f (x)dx b − an

n−1i=0

f xi + xi+1

2

poart˘ a numele de formula dreptunghiurilor, denumire care provine din interpretarea geometric˘ a pentru f ≥ 0. ˆ Intr-adev˘ ar, dac˘ a Di este drept-

unghiul care are ca baz˘ a intervalul [xi, xi+1] si ın alt ime f

xi + xi+1

2

,

atunci

b a

f (x)dx (adic˘ a aria ”de sub graficul” lui f ) poate fi aproximat˘ a cu

suma tuturor ariilor dreptunghiurilor Di, i = 0,...,n − 1 (vezi figura 5.1).

Exemplul 14 S˘ a se calculeze valoarea aproximativ˘ a a integralei

1 0

dx

1 + x2

cu o eroare ≤ 1

40, utilizˆ and formula dreptunghiurilor.

Fie f : [0, 1] → IR, f (x) =1

1 + x2. Atunci f ∈ C 1([a, b]), iar f (x) =

−2x

(1 + x2)2. Vom ar˘ ata acum c˘ a f ∞ ≤ 1. Avem succesiv

2x

(1 + x2)2≤ 1 ⇔

2x ≤ (1 + x2)2, ∀ x ∈ [0, 1]. Fie acum funct ia ϕ(x) = x4 + 2x2 − 2x + 1 ≥0, ∀ x ∈ [0, 1], de unde ϕ(x) = 2(2x3 + 2x − 1) si folosind sirul lui

Rolle se obt ine c˘ a ϕ admite o singur˘ a r˘ ad˘ acin˘ a real˘ a ın (0,1

2 ) pentru c˘ a ϕ(0)ϕ( 1

2 ) < 0 si ϕ este continu˘ a. Deci

ϕ(x) > ϕ

1

2

=

1

2> 0, ∀ x ∈

1

2, 1

,

iar pentru

x ∈

0,1

2

, 2x ≤ 1 ≤ (1 + x2)2 ⇒ ϕ(x) ≥ 0, ∀ x ∈

0,

1

2

si deci

ϕ(x)≥

0,∀

x∈

[0, 1]⇒

f ∞ ≤

1.

Determin˘ am acum n astfel ınc at (1 − 0)2

4n1 ≤ 1

40si obt inem n ≥ 10. Lu˘ am

n = 10 ın (5.6) si obt inem 1

0f (x)dx 1

10

f

1

20

+ f

3

20

+ ... + f

19

20

= 0.78560636

49



cu o eroare mai mic˘ a decˆ at 1

40. Cum ıns˘ a

1

0

dx

1 + x2= arctg(x)/1

0 =π

4

rezult˘ a c˘ a ıl putem aproxima pe π prin 4 · 0.78560636 = 3, 14242534 cu o

eroare mai mic˘ a decˆ at 1

40dat˘ a de metoda dreptunghiurilor.

Teorema 22 (Formula trapezelor) Fie f : [a, b] → IR, f ∈ C 2([a, b])si ∆ : (a = x0 < ... < xn = b) o diviziune echidistant˘ a pentru [a, b]. Definim num˘ arul σ2(f, n) prin

σ2(f, n) =b − a

n

n−1i=0

f (xi) + f (xi+1)

2. (5.9)

Atunci σ2(f, n) aproximeaz˘ a ba

f (x)dx cu o eroare majorat˘ a de

(b − a)3

12n2f ∞.

Observatia 39 Formula (5.9) este exact˘ a pentru f ∞ = 0, adic˘ a pentru funct iile polinomiale de grad cel mult 1.

Figura 5.2: Metoda trapezelor

Observatia 40 Relat ia ba

f (x)dx b − a

n

ni=0

f (xi) + f (xi+1)

2

50



poart˘ a numele de formula trapezelor, denumire care provine de asemenea

din interpretarea geometric˘ a pentru f ≥

0. b

af (x)dx se aproximeaz˘ a cu

suma ariilor tuturor trapezelor T i, cu A(T i) =b − a

n· f (xi) + f (xi+1)

2(vezi

figura 5.2).

Exemplul 15 S˘ a se aproximeze integrala

1

0

dx

1 + xcu o eroare mai mic˘ a

decˆ at 1

600, folosind formula trapezelor.

Fie f : [0, 1] → IR, f (x) =1

1 + x. Atunci | f ”(x) |= 2

(1 + x)3≤ 2, ∀ x ∈

(0, 1]. Determin˘ am apoi pe n astfel ınc at (1 − 0)3

12n2 ·2

≤1

600si obt inem

n ≥ 10. Astfel 1

0f (x)dx 1

10

f (0) + f (1)

2+ f

1

10

+ ... + f

9

10

= 0.69377136,

deci 1

0

dx

1 + x≈ 0.69377136

cu o eroare mai mic˘ a decˆ at ≤ 1

100. Cum ıns a

1

0

dx

1 + x= ln2 rezult˘ a c˘ a

ln 2 poate fi aproximat prin valoarea 0.69377136 cu o eroare mai mic˘ a decˆ at

1600

dat˘ a de metoda trapezelor.

Teorema 23 (Formula Simpson) Fie f : [a, b] → IR, f ∈ C 4([a, b]), ∆ :(a = x0 < .. < xn = b) o diviziune echidistant˘ a pe [a, b] si num arul σ3(f, n)definit prin

σ3(f, n) =b − a

6n

n−1i=0

f (xi) + 4f

xi + xi+1

2

+ f (xi+1)

. (5.10)

Atunci σ3(f, n) aproximeaz˘ a pe b

a

f (x) dx cu o eroare majorat˘ a de

(b − a)5

2880n4f (1)∞. (5.11)

Observatia 41 Aceast˘ a formul˘ a este exact˘ a pentru f (4)∞ = 0 adic˘ a pen-tru funct iile polinomiale de grad cel mult 3.

51



Figura 5.3: Metoda Simpson

Observatia 42 Relat ia ba

f (x) dx b − a

6n

n−1i=0

f (xi) + 4f

xi + xi+1

2

+ f (xi+1)

(5.12)

poart˘ a denumirea de formula Simpson sau a parabolelor si provine din interpretarea geometric˘ a pentru f ≥ 0 (vezi figura 5.3). Aria trapezului curbiliniu ABCD este aproximat˘ a cu aria trapezului curbiliniu ABCFDce are drept una din laturile neparalele parabola CF D (de ecuat ¸ie y =

ax2 + bx + c) care este unic˘ a cu proprietatea c˘ a β

α(ax2 + bx + c) dx =

β − α

6

f (α) + 4f

α + β

2

+ f (β )

. (5.13)

Aici, [α, β ] = [xi, xi+1].

Exemplul 16 Utilizˆ and formula Simpson, s˘ a se aproximeze integrala 1

0

dx

x + 1, cu o eroare mai mic˘ a decˆ at

1

10000.

Fie f : [0, 1] → IR, f (x) =1

x + 1cu f (4)(x) =

4!

(1 + x)5. Atunci f (4)∞ =

supx∈[0,1]

|4!(1+x)−5| = 4! = 24. Determin˘ am apoi n astfel ıncˆ at (1 − 0)5

2880n4 ·24 <

1

104si obt inem n ≥ 4. Atunci 1

0f (x) dx 1 − 0

6 · 4

f (0) + f (1) + 2

f

1

4

+ f

2

4

+ f

3

4

+

52



+4

f

1

8

+ f

3

8

+ f

5

8

+ f

7

8

= 0.69315449

deci 10 dx1 + x

0.69315449 cu o eroare majorat˘ a de 10−4. Cum ıns˘ a 1

0

dx

1 + x= ln2 rezult˘ a c˘ a ln 2 poate fi aproximat prin valoarea 0.69315459

cu o eroare mai mic˘ a decˆ at 10−4 dat˘ a de metoda Simpson.

Observatia 43 Metoda Simpson se poate aplica ın urm˘ atoarea form˘ a.

1. Se alege n ∈ IN ∗ par;

2. Se ımparte intervalul [a, b] ın n subintervale de lungime egal˘ a. Pentru

h =b − a

n, xi = a + ih, i = 0, . . . , n avem

ba

f (x) dx =

n/2i=1

x2ix2i−1

f (x) dx n/2i=1

h

3[f (x2i−2) + 4f (x2i−1) + f (x2i)] ,

∀ f ∈ C 4([a, b]).

Observatia 44 1. Formula anterioar˘ a este folosit˘ a ın aplicat ii sub forma

ba

f (x) dx h

3

f (a) + 2

(n/2)−1i=1

f (x2i−1) + 4

n/2i=1

f (x2i−1) + f (b)

,

(5.14)care se numeste formula Simpson compozita.

2. Ca si ın cazul formulei Simpson, se poate arat˘ a (vezi [ 3 ], [ 16 ]) c˘ a eroarea ın cazul formulei Simpson compozite este

b − a

180h4f (4)∞.

3. ˆ In acelasi mod se pot obt ine formula trapezelor compozita, res-pectiv formula dreptunghiurilor compozita.

5.1.2 Aproximarea integralelor pe IR2

In sectiunea anterioara am prezentat formule de tip Newton-Cotes pentruaproximarea integralelor pe IR. Formule analoage exista si pentru aproxi-marea integralelor duble. Astfel, pentru a aproxima integrala dubla

D

f (x, y) dx dy

53



unde f : D −→ IR este (cel putin) continua si

D = (x, y) ∈ IR2

/a ≤ x ≤ b si c ≤ y ≤ d.

se poate folosi orice formula descrisa anterior. Pentru a aproxima prinformula Simpson compozita, se aleg numerele pare n si m si se consideradiscretizarea x0 = a < x1 < ... < xn = b, respectiv y0 = c < y1 < ... <

yn = d, unde h =b − a

n, iar k =

d − c

m. Conform teoremei lui Fubini (vezi

[15], [26]) avem D

f (x, y) dx dy =

ba

dc

f (x, y)dy

dx.

Atunci, aproximam mai ıntai dc f (x, y) dy. Astfel, din relatia (5.14) avem dc

f (x, y)dy =k

3

f (x, y0) + 2

(m/2)−1 j=1

f (x, y2 j) + 4

m/2 j=1

f (x, y2 j−1) + f (x, ym)

deci b

a

dc

f (x, y)dy dx =k

3

ba

f (x, y0) dx + 2

(m/2)−1 j=1

ba

f (x, y2 j) dx+

+4

m/2 j=1

ba

f (x, y2 j−1) dx + ba

f (x, ym) dx .

Acum se aproximeaza fiecare integrala ın parte, luand xi = a + ih, i =0, . . . , n. Pentru orice j = 0, . . . , m avem

ba

f (x, y j)dy =h

3

f (x0, y j) + 2

(n/2)−1i=1

f (x2i, y j) + 4

n/2i=1

f (x2i−1, y j)+

+f (xn, y j)] , deci

ba

dc

f (x, y) dy dx kh3f (x0, y0) + 2

(n/2)−

1i=1

f (x2i, y0)+

+4

n/2i=1

f (x2i−1, y0) + f (xn, y0)

+ 2

(m/2)−1 j=1

f (x0, y2 j)+

54



+2

(m/2)−1

j=1

(n/2)−1

i=1

f (x2i, y2 j) + 4

(m/2)−1

j=1

n/2

i=1

f (x2i−1, y2 j)+

+

(m/2)−1 j=1

f (xn, y2 j)

+ 4

m/2 j=1

f (x0, y2 j−1) + 2

m/2 j=1

(n/2)−1i=1

f (x2i, y2 j−1)+

+4

m/2 j=1

n/2i=1

f (x2i−1, y2 j−1) +

m/2 j=1

f (xn, y2 j−1)

+ [f (x0, ym)+

+2

(n/2)−1i=1

f (x2i, ym) + 4

n/2i=1

f (x2i−1, ym) + f (xn, ym)

.

Observatia 45 ˆ In formula anterioar˘ a, termenii din sum˘ a se pot grupa si aceasta se poate pune sub forma

hk

9

ni=0

m j=0

wi,jf (xi, y j).

Exemplul 17 Aproximat i I =2.0

1.4

1.5 1.0

ln(x + 2y) dy dx, folosind formula

Simpson compozit˘ a pentru n = 4 si m = 2.

Avem h =0.6

4= 0.15, k =

0.5

2= 0.25, deci

I 0.15 · 0.25

9·

4i=0

2 j=0

wi,j ln(xi + 2y j) = 0.4295524387

; valoarea exact˘ a (cu 10 zecimale) a integralei este 0.4295546265, ceea ce ınseamn a c˘ a eroarea este de ordinul 10−6.

Observatia 46 ˆ In mod analog pot fi deduse funct iile de cuadratur˘ a com-pozite pentru metoda trapezelor pentru aproximarea integralelor duble si funct iile de cuadratur˘ a compozite (trapeze si Simpson) pentru aproximarea integralelor multiple (n ≥ 3) (pentru detalii vezi lucr˘ arile [ 3 ] , [ 16 ] , [ 27 ] ).

5.2 Formule de cuadratura cu noduri Gauss

5.2.1 Formule de tip Cebısev

Principalul dezavantaj al formulelor de integrare de tip Newton-Cotes esteacela ca gradul polinomului de interpolare folosit creste odata cu numarul

55



de puncte din diviziune si ın plus, punctele sunt echidistante, deci, implicita priori fixate, fara legatura cu unele proprietati functiei de sub integrala.

Aceste inconveniente pot fi eliminate prin considerarea formulelelor detip Cebısev. In acest sens consideram urmatoarea problema:pentru f : [−1, 1] → IR continua, sa se determine punctele t1, ..., tn ∈[−1, 1] astfel ıncat formula 1

−1f (t)dt = B

ni=1

f (ti) = σ4(f, n) (5.15)

sa fie adevarata pentru orice polinom de grad cel mult n. Scriind (5.15)pentru polinoamele din baza canonica 1, x , x2,..,xn, pentru f = 1 rezulta

B =2

n, iar

σ4(f, n) = 2n

ni=1

f (ti), (5.16)

unde t1,..,tn sunt solutiile (reale) ale sistemului de ecuatii neliniare

t1 + t2 + ... + tn = 0

tk1 + tk

2 + ... + tkn =

n

2

1 − (−1)k+1

k + 1.....................................

tn1 + tn

2 + .. + tnn =

n

2(n + 1)

1 − (−1)n+1

(5.17)

Observatia 47 Punctele t1, ...,tn nu depind de funct ia f .Observatia 48 S. Bernstein a demonstrat (vezi [ 6 ]) c˘ a pentru n = 8 si n ≥ 10, sistemul (5.17) are si solut ii complexe, deci singurele cazuri posibile sunt

n = 2, 3, 4, 5, 6, 7, 9. (5.18)

Aceste valori sunt deja calculate si se g˘ asesc ın tabele. De exemplu, pentru n = 2 r˘ ad˘ acinile sunt x1 = 0.5773502692, x2 = −0.5773502692, iar pentru n = 3, x1 = 0.7745966692, x2 = 0.0000000000, x3 = −0.7745966692.

Observatia 49 Formulele de tip Cebısev sunt destul de precise si sunt utile cˆ and avem de calculat foarte multe integrale pe acelasi domeniu (vezi de ex.

discretiz˘ arile de tip element finit pentru ecuat ii cu derivate part iale, [ 4 ]).

Observatia 50 Extinderea considerat iilor anterioare de la [−1, 1] la un interval arbitrar [a, b] ⊂ IR, se face cu ajutorul schimb˘ arii de variabil˘ a

x =b + a

2+

b − a

2t. (5.19)

56



5.2.2 Formule bazate pe polinoame Legendre

O alta metoda de tip Gauss se bazeaza pe utilizarea radacinilor poli-noamelor Legendre P 0(x), P 1(x),...,P n(x),..., caracterizate de urmatoare-le proprietati

1. oricare ar fi n, P n este polinom de gradul n.

2.

1

−1P (x) · P n(x) dx = 0 pentru orice polinom P (x) de gradul cel mult

n − 1.

3. Radacinile polinoamelor P n, n ≥ 0 sunt distincte, sunt ın intervalul(−1, 1) si simetrice fata de origine.

Observatia 51 Primele cˆ ateva polinoame Legendre sunt

P 0(x) = 1, P 1(x) = x, P 2(x) = x2 − 1

3, P 3(x) = x3 − 3

5x

P 4(x) = x4 − 6

7x2 +

3

35.

Este adevarat urmatorul rezultat (pentru demonstratie vezi [3]).

Teorema 24 Dac˘ a x1, x2, ..., xn sunt r˘ ad˘ acinile polinomului Legendre P n(x) si dac a pentru orice i = 1, . . . , n definim numerele

ai = 1

−1

n j=1, j=i

x−

x j

xi − x jdx,

atunci pentru orice polinom de grad maxim 2n − 1 avem 1

−1P (x) dx =

ni=1

aiP (xi).

5.2.3 Formule de tip Gauss pentru integrale duble

Ca si ın cazul integrarii cu noduri Newton-Cotes se pot folosi noduri detip Gauss si pentru aproximarea integralelor multiple (n

≥2). Ca noduri

Gauss se folosesc ın acest caz radacinile polinomului Legendre prezentatanterior.

Exemplul 18 Pentru a aproxima integrala din exemplul (17) adic˘ a

I =

2.0

1.4

1.5

1.0ln(x + 2y) dy dx

57



transform˘ am mai ınt ai dreptunghiul D = [1.4, 2.0] × [1.0, 1.5] ın dreptunghiul D =

(u, v)/

−1

≤u

≤1 si

−1

≤v

≤1

; astfel rezult˘ a

(vezi si (5.19))

u =1

2 − 1.4(2x − 1.4 − 2) si v =

1

1.5 − 1(2y − 1 − 1.5),

de unde obt inem x = 0.3u + 1.7, y = 0.25v + 1.25 si

I = 0.3 · 0.25

1

−1

1

−1ln(0.3u + 0.5v + 4.2)dvdu.

Pentru integrarea cu noduri Gauss pentru n = 3 vom folosi

u0 = v0 = x3 = −0.7745966692

u1 = v1 = x2 = 0

u2 = v2 = x1 = −x3 = 0.7745966692,

deci

I = 0.0753

i=1

3i=1

a3,ia3,j ln(0.3xi + 0.5x j + 4.2) = 0.4295545313,

care aproximeaz˘ a integrala dat˘ a cu o eroare de ordinul 10−9, spre deosebire de metoda Simpson compozit˘ a, pentru care aveam doar 10−6.

Observatia 52 Dac˘ a ın aceast a sect iune si ın 5.1.2 ne-am referit la aproxi-marea integralelor pe domenii dreptunghiulare, exist˘ a cazuri cˆ and ne intere-seaz˘ a aproximarea integralelor pe triunghiuri. Acestea apar, de exemplu, la discretizarea ecuat iilor cu derivate part iale prin metoda elementului finit (vezi [ 4 ]).

Astfel, fie triunghiul T cu varfurile Ai(xi, yi), i = 1, 2, 3 si fie M, N ,respectiv P mijlocul laturii A1A2, A1A3, respectiv A2A3, iar G baricentrul(centrul de greutate) al tringhiului. Se stie ca

xM =x1 + x2

2, yM =

y1 + y2

2,

xN = x1 + x32 , yN = y1 + y3

2 ,

xP =x2 + x3

2, yM =

y2 + y3

2

si

xG =x1 + x2 + x3

3, yG =

y1 + y2 + y3

3.

58



In [Ciarlet] sunt prezentate mai multe formule de cuadratura care sunte- xacte pentru anumite clase de functii. Daca notam cu A(T ) aria tri-

unghiului T si cu f functia de integrat pe triunghiul T , avem la dispozitieurmatoarele formule de cuadraturaFormula 1.

T

f (x, y) dx dy A(T )

3[f (A1) + f (A2) + f (A3)] = A(T ) · f (G). (5.20)

Formula este exacta pentru polinoame de grad cel mult 1.

Exemplul 19 Dac˘ a f (x, y) = 1, ∀ x,y, atunci T

f (x, y) dx dy (care

reprezint˘ a, de fapt aria A(T ) a triunghiului T ), este aproximat˘ a prin (5.20)

de valoarea A(T )

3 (1 + 1 + 1) = A(T ) · 1 = A(T ). Se confirm˘ a astfel exacti-tatea formulei (5.20) ment ionat˘ a anterior.

Formula 2. T

f (x, y) dxdy A(T )

3[f (M ) + f (N ) + f (P )] . (5.21)

Formula este exacta pentru polinoame de grad cel mult 2.

Exemplul 20 Dac˘ a f (x, y) = xy, ∀ x,y, iar A1(0, 0), A2(0, 1), A3(0, 0),

atunci T f (x, y) dxdy ( a c˘ arei valoare exact˘ a este dat˘ a de

1

0 1−x

0 xy dy dx =

1/24), se aproximeaz˘ a prin (5.21) cu valoarea A(T )

3

0 · 1

2+

1

2· 0 +

1

2· 1

2

=

A(T )/12 = 1/24. Se confirm˘ a astfel exactitatea formulei (5.21) pentru poli-noame de grad cel mult 2.

Formula 3. T

f (x, y) dx dy A(T )

60[3 · (f (A1) + f (A2) + f (A3)) +

8 · (f (M ) + f (N ) + f (P )) + 27 · f (G)]. (5.22)Formula este exacta pentru polinoame de grad cel mult 3.

Exemplul 21 Cu triunghiul T ales ca ın exemplele anterioare si funct ia f (x, y) = x2y, se constat˘ a c˘ a formula (5.22) este de asemenea exact˘ a pentru T

f (x, y) dx dy.

59



Vom remarca faptul ca ın formulele prezentate mai sus, cu cat preciziaeste mai mare, cu atat numarul de puncte ın care trebuie calculata va-

loarea functiei de integrat este mai mare. De aceea, ın aplicat ii la metodaelementului finit (unde numarul de triunghiuri din discretizare este foartemare: 105 − 106, vezi [4]), trebuie gasita calea optima ın ceea ce privesteefortul de calcul.

5.3 Exercitii

1. Sa se aproximeze

1

0

x

x2 + 1dx prin metoda dreptunghiurilor cu o

eroare mai mica sau egala decat1

400.

2. Sa se aproximeze 1

0ln(x2 + 1) dx prin metoda trapezelor cu o eroare

mai mica sau egala decat1

2400.

3. Sa se aproximeze

1/2

0

2x − 1

x2 + 4dx prin metoda dreptunghiurilor cu o

eroare mai mica sau egala decat1

800.

4. Sa se aproximeze 1/2

0

sin x

x

dx cu o eroare mai mica a decat 10−6

(a) dezvoltand mai ıntai functia continua f :

0,

1

2

→ IR,

f (x) =

sin x

x, x = 0

1 , x = 0ın serie Taylor si apoi integrand

termen cu termen (seria obtinuta fiind uniform convergenta pe0,

1

2

);

(b) prin metoda dreptunghiurilor;

(c) prin metoda trapezelor;

(d) prin metoda Simpson.

5. Sa se aproximeze

3

2

4x + 1

(x + 1)3dx prin metoda trapezelor cu o eroare

mai mica sau egala decat1

2700.

60



6. Sa se aproximeze

3

2

x2 + x + 2

(x − 1)2(x + 1)2dx prin metoda dreptun-

ghiurilor cu o eroare mai mica sau egala decat 13000

.

7. Aproximati integralele urmatoare folosind metoda trapezelor

(a)

1.5

1x2 ln x dx

(b)

0.35

0

2

x2 − 4dx

(c) π/4

0 x sin x dx

pentru n = 2.

8. Reluati exercitiul 7 folosind metoda Simpson cu n = 3.

9. Aproximati urmatoarea integrala folosind noduri Gauss pentru n = 2si apoi comparati valorile obtinute cu valoarea exacta

(a)

1

0x2e−x dx,

(b)

π/4

0(cos x)2 dx,

(c) π/4

0

e3x sin2x dx.

10. Reluati exercitiul 9 pentru n = 3.

11. Pentru m = n = 4, aproximati integralele urmatoare folosind formulaSimpson compozita:

(a)

2.5

2.1xy2 dy dx;

(b)

0.5

0ey−x dy dx.

12. Reluati exercitiul 11 pentru n = m = 2 folosind noduri Gauss.

61



Capitolul 6

Aproximarea solutiilor

ecuatiilor diferentiale

ordinare

Fie problema cu valori initialey = f (t, y)y(a) = α,

(6.1)

unde f : [a, b] × I → IR este continua, iar I ⊂ IR este un interval. Prinsolutie a problemei (6.1) se ıntelege o aplicatie y : [a, b] → I de clasa C 1 ceverifica relatiile din (6.1),

∀t

∈[a, b].

Definitia 10 Aplicat ia f se numeste Lipschitziana ın raport cu a doua variabil˘ a dac˘ a exist˘ a constanta L ≥ 0 astfel ınc at

|f (t, y1) − f (t, y2)| ≤ L |y1 − y2| , ∀ (t, yi) ∈ [a, b] × I, i = 1, 2. (6.2)

Observatia 53 Din teoria general˘ a a a ecuat ¸iilor diferent ¸iale (vezi [ 3 ] ,[ 22 ] si referintele lor) se stie c˘ a o condit ie suficient˘ a pentru a avea loc (6.2)este ca

∂f

∂y(t, y)

≤ L, ∀ (t, y) ∈ [a, b] × I. (6.3)

Teorema 25 Fie f : [a, b]×

IR→

IR continu˘ a si Lipschitzian˘ a ın raport cu a doua variabil˘ a. Atunci problema (6.1) are solut ie unic˘ a.

Definitia 11 Problema (6.1) se numeste corect pusa dac˘ a are solut ie,aceasta este unic˘ a si continu a ın raport cu datele init iale, deci este stabila la perturbat ii (adic˘ a, perturbat ii mici ale datelor init iale, f si α determin˘ a perturbat ii mici ale solut iei).

62



Teorema 26 ˆ In condit iile Teoremei 25, problema (6.1) este corect pus˘ a.

Observatia 54 Necesitatea introducerii not iunii de problem˘ a corect pus˘ a se datoreaz˘ a faptului c˘ a, ın mare parte, problemele de tipul (6.1) mode-leaz˘ a niste fenomene fizice reale. Astfel, expresiile pentru f si α cont in unele date/coeficient i supuse erorilor de modelare (m˘ asur˘ atori, aproxim˘ ari de modelare, etc).

Observatia 55 Cu except ia unor cazuri particulare, problemele de tipul (6.1) nu pot fi rezolvate exact (adic˘ a s˘ a g˘ asim o expresie analitic˘ a pentru solut ie). De aceea se apeleaz˘ a de rezolvarea lor numeric˘ a, adic˘ a obt inerea unei aproxim˘ ari a solut iei exacte.

6.1 Metode de tip Taylor

Vom prezenta ın continuare o clasa generala de metode pentru aproximareasolutiei problemei (6.1). Fie ∆ = (a = t0 < t1 < ... < tN = b) o diviziuneechidistanta a lui [a, b], adica

h =b − a

N , ti = a + ih, i = 0, . . . , N . (6.4)

Pentru h0 > 0 fixat, definim D = [a, b] × IR × [0, h0] ⊆ IR3 si consideram ofunctie φ : D → IR continua si Lipschitziana ın raport cu a doua variabila,adica ∃ L > 0 astfel ıncat ∀ (t, wi, h) ∈ D, i = 1, 2

|φ(t, w1, h) − φ(t, w2, h)| ≤ L |w1 − w2| . (6.5)

Presupunem ın plus ca φ satisface conditia de consistenta

φ(t, y(t), 0) = f (t, y(t)), ∀ t ∈ [a, b], (6.6)

unde y : [a, b] → IR este solutie exacta pentru problema (6.1). Consideramurmatoarea metoda generala de aproximare a solutiei problemei (6.1)

w0 = αwi+1 = wi + h · φ(ti, wi, h), i = 1, . . . , N .

(6.7)

Pentru metoda (6.7) definim

E i = y(ti) − wi = yi − wi, E = maxi=1,N

|E i| (6.8)

care se numesc eroare locala absoluta, respectiv eroare globala ab-soluta si masoara diferenta dintre aproximatiile wi si valorile exacte y(ti),

63



oferind o informatie clara despre eficienta aproximarii date de (6.7). Con-sideram de asemenea valorile

τ i =y(ti) − y(ti−1)

h− φ(ti−1, y(ti−1), h), i = 1, . . . , N , (6.9)

siτ = max

i=1,N |τ i| (6.10)

numite eroarea locala de trunchiere, respectiv eroarea globala detrunchiere. Acestea masoara abaterea de la egalitatea ce s-ar obtine prinintroducerea valorilor exacte y(ti) ın (6.7).

Definitia 12 Spunem c˘ a metoda (6.7) este de ordin λ dac˘ a exist˘ a c0 si λ constante pozitive, independente de N astfel ınc at

E ≤ c0 · hλ. (6.11)

Observatia 56 Pentru anumite clase de metode de tip (6.7) este relativ simplu s˘ a g˘ asim constantele c, β independente de N astfel ınc at

τ ≤ c · hβ . (6.12)

Urmatorul rezultat ne asigura ca pentru a avea o conditie de tipul (6.11) esuficient sa fie ındeplinita o conditie de tipul (6.12).

Teorema 27 Dac˘ a are loc (6.12), atunci avem

E ≤ hβ ·c

L eL(b−a) − 1 , (6.13)

adic˘ a metoda (6.7) este de ordin β .

Demonstratie. Din (6.7) obtinem

y(ti+1) − wi+1 = yi+1 − wi+1 = yi+1 − wi − h · φ(ti, wi, h).

Cum ınsa din (6.11) avem

yi+1 = h · τ i+1 + h · φ(ti, wi, h),

rezulta

yi+1 − wi+1 = yi − wi + h · τ i+1 + h [φ(ti, yi, h) − φ(ti, wi, h)] .

Fie acum ai = |E i| = |yi − wi| , i = 0, . . . , N . Folosind si relatia (6.12)rezulta

ai+1 ≤ ai + hβ +1 · c + h · L · ai = ai(1 + hL) + c · hβ +1,

64



Teorema 28 Dac˘ a M ≥ 0 este astfel ıncˆ at

y(n+1)

(t) ≤ M, ∀ t ∈ [a, b] (6.19)

si dac a f, f ,..,f (n−1) sunt Lipschitziene ın raport cu a doua variabil˘ a f (k)(t, y1) − f (k)(t, y2) ≤ λ · |y1 − y2| , ∀ (t, yi) ∈ [a, b] × IR, (6.20)

i = 1, 2, ∀ k = 0, . . . , n − 1 atunci, pentru metoda Taylor (6.18)avem evaluarea

E ≤ hn · M

n! · λ

e2λ(b−a) − 1

. (6.21)

Observatia 57 Relat ia (6.21) ne spune c˘ a metoda Taylor este de ordin n.Pentru n = 1 obt inem metoda Euler

w0 = α0

wi+1 = wi + h · f (ti, wi), i = 1, . . . , N − 1.(6.22)

care este de ordin 1.

Evaluarea (6.21) ne spune ca, pentru valori mari ale lui n, metoda Tayloreste foarte eficenta. Dar, practic aceasta eficienta poate fi estompata sauchiar anulata de numarul mare si complexitatea calculelor necesare evaluariicantitatii T n(ti, wi) din (6.17).

6.2 Metode de tip Runge-KuttaPentru a elimina inconvenientul precizat ın sectiunea anterioara (cu pastrareaunui ordin suficient de bun) au fost create metodele de tip Runge-Kutta.S-a pus problema (de exemplu) sa construim o metoda de tipul (6.7) caresa aiba un ordin mai mare ca 1 (metoda Euler), dar fara a folosi ın calculullui wi derivatele functiei f . In continuare vom deduce o astfel de metodapentru n = 2. Din (6.17) rezulta

T 2(t, y) = f (t, y) +h

2!f (t, y)

sau utilizand (6.15)

T 2(t, y) = f (t, y) +h

2· ∂f

∂t(t, y) +

h

2· ∂f

∂y(t, y) · f (t, y). (6.23)

Vom aproxima T 2(t, y) printr-o expresia de forma

A · f (t + α, y + β ), A, α, β ∈ IR (6.24)

66



cu o eroare de ordinul O(h2). Dezvoltand (6.24) ın serie Taylor de ordinul2 obtinem

A · f (t + α, y + β ) = A · f (t, y) + Aα · ∂f

∂t(t, y)+

Aβ · ∂f

∂y(t, y) + A · R(t + α, y + β ), (6.25)

unde

R(t + α, y + β ) =α2

2· ∂ 2f

∂t2(ξ, η) + α · β · ∂ 2f

∂t∂y(ξ, η) +

β 2

2· ∂ 2f

∂y2(ξ, η) (6.26)

cu ξ, η ıntre t si t + α, y si y + β , respectiv. Identificand coeficientiiderivatelor partiale ın (6.23) si (6.25) rezulta (ın mod unic) valorile A = 1,α = h

2 , β = h2 · f (t, y). Se obtine astfel metoda punctului de mijloc

w0 = α

wi+1 = wi + h · f

ti +

h

2, wi +

h

2· f (ti, wi)

,

(6.27)

i = 1, . . . , N − 1. Alte metode de ordinul 2 se pot obtine daca ın (6.24) sefolosesc combinatii de forma

A · f (t, y) + B · f (t + α, y + β · f (t, y)). (6.28)

Rezulta astfel metoda Euler modificataw0 = α

wi+1 = wi +h

2· [f (ti, wi) + f (ti + h, wi + h · f (ti, wi))] .

(6.29)

sau metoda lui Heun.

w0 = α

wi+1 = wi + h4

· f (ti, wi) + 3f (ti + 23

h, wi + 23

h · f (ti, wi)) ,

i = 1, . . . , N − 1. Daca se folosesc expresii mai complicate ın (6.28) se obtinmetode de tip Runge-Kutta de ordin mai mare. Varianta clasica a acesteimetode (care este de ordinul 4 daca y ∈ C 5([a, b])) este urmatoarea (pentru

67



detalii si dezvoltari vezi [???])

w0 = αk1 = h · f (ti, wi)

k2 = h · f (ti +h

2, wi +

k1

2) (6.30)

k3 = h · f (ti +h

2, wi +

k2

2)

k4 = h · f (ti + h, wi + k3)

wi+1 = wi +k1 + 2k2 + 2k3 + k4

6, i = 1, . . . , N − 1.

Observatia 58 Metodele prezentate ın acest capitol pot fi extinse si la sis-teme de ecuat ii diferent iale ordinare de tipul (6.1) sau chiar la ecuat ii si

sisteme de ecuat ii de ordin superior (vezi [ 3 ]).

6.3 Exercitii

1. Folosind metoda lui Euler, aproximati solutia problemei cu valorileinitiale urmatoare

(a) y = te3t − 2y, t ∈ [0, 1], y(0) = 0 si h = 0.5;

(b) y = cos 2t + sin 3t, t ∈ [0, 1], y(0) = 1 si h = 0.25;

(c) y = 1 +y

t+

y2

t2, t ∈ [1, 3], y(1) = 0 si h = 0.2.

2. Calculati eroarea dintre solutia exacta a problemei de la exercitiul 1si cea data de metoda lui Euler la fiecare pas.

3. Folosind metoda lui Taylor de ordinul 2 aproximati solutia problemeicu valoarile initiale

(a) y = te3t − 2y, t ∈ [0, 1], y(0) = 0 si h = 0.5;

(b) y = 1 + (t − y)2, t ∈ [2, 3], y(2) = 1 si h = 0.5;

(c) y = cos 2t + sin 3t, t ∈ [0, 1], y(0) = 1 si h = 0.5.

4. Folosind metoda lui Euler modificata aproximati

(a) y = te3t − 2y, t ∈ [0, 1], y(0) = 0 si h = 0.5

(b) y = 1 + (t − y2), t ∈ [2, 3], y(2) = 1 si h = 0.5

(c) y = cos 2t + sin 3t, t ∈ [0, 1], y(0) = 1si h = 0.25

si comparati cu valoarea exacta.

68



5. Sa se aproximeze solutia urmatorului sistem de ecuatii diferentialeordinare folosind metoda lui Runge-Kutta; comparati rezultatul obtinut

cu solutia clasica.

(a)

u1 = 3u1 + 2u2 − (2t2 + 1)e2t , u1(0) = 1u2 = 4u1 + u2 + (t2 + 2t − 4)e2t , u2(0) = 1

pentru h = 0.2 si t ∈ [0, 1]

(b)

u1 = u2 − u3 + t , u1(0) = 1u2 = 3t2 , u2(0) = 1u3 = u2 + e−t , u3(0) = −1

pentru h = 0.1 si t ∈ [0, 1]

6. Folosind metoda Runge-Kutta pentru sisteme, aproximati solutia urmatoarelorecuatii diferentiale de ordin superior:

(a)

y − 2y + y = tet − t , t ∈ [0, 1]y(0) = y(0) = 0

si h = 0.1

(b)

y + 2y − y − 2y = et , t ∈ [0, 3]y(0) = 1, y(0) = 2, y(0) = 0

cu h = 0.2

69



Capitolul 7

Metode directe pentru

sisteme liniare

Se considera sistemul liniar

Ax = b. (7.1)

cu A ∈ Mm×n, b ∈ IRm si x ∈ IRn. Pentru a rezolva sistemul (7.1) existadoua clase mari de metode:

1. metode directe, ın care solutia sistemului se obtine dupa un numarfinit de pasi, posibil de evaluat a priori;

2. metode iterative, ın care solutia este limita unui sir de aproximatii.

Metodele iterative vor fi studiate ın Capitolul 8. In prezentul capitol nevom ocupa de prezentarea unor metode directe de rezolvare.

7.1 Metoda eliminarii a lui Gauss

Inainte de a prezenta metoda eliminarii a lui Gauss (pe scurt, EG), vomintroduce unele notiuni si notatii.

Definitia 13 O matrice L ∈ Mn se numeste inferior triunghiulara (pe scurt, inf.∆) dac˘ a

lij = 0, ∀ j > i, i = 1, . . . , n . (7.2)

Definitia 14 O matrice U ∈ Mn se numeste superior triunghiulara(pe scurt, sup.∆) dac˘ a

uij = 0, ∀ j < i, i = 1, . . . , n . (7.3)

70



Definitia 15 O matrice D ∈ Mn se numeste diagonala (pe scurt, diag)dac˘ a

dij = 0, ∀ j = i, i = 1, . . . , n . (7.4)

Vom nota cu Li linia i din sistemul (7.1), adica

ai1x1 + ai2x2 + · · · + ainxn = bi. (7.5)

Definitia 16 Se numeste transformare elementara a sistemului (7.1)oricare din urm˘ atoarele operat ii

• linia i adunat˘ a cu linia j se pune ın locul liniei i

• linia i ınmult it˘ a cu scalarul λ se pune ın locul liniei i• liniile i si j comut˘ a ıntre ele

Vom nota aceste transformari prin

Li + L j −→ Li; λ · Li −→ Li; Li ←→ L j. (7.6)

Avem urmatoarele rezultate (pentru demonstratii vezi [23], [24]).

Propozitia 2 (i) Produsul a dou˘ a matrice inf.∆ (respectiv sup.∆) este omatrice inf.∆ (respectiv sup.∆).(ii) Dac˘ a o matrice inf.∆ (respectiv sup.∆) este inversabil˘ a, atunci inversa

sa este, de asemenea, o matrice inf.∆ (respectiv o matrice sup.∆).(iii) Mult imea solut iilor sistemului (7.1) este invariant˘ a la transform˘ arile elementare (7.6).

Vom prezenta ın continuare metoda EG. Aceasta consta ın urmatorii pasi:

1. Transformarea sistemului Ax = b ıntr-un sistem Ax = b, cu A matricesuperior triunghiulara si inversabila (atunci cand este posibil), adica

a11x1 + . . . + a1nxn = b1

a21x1 + . . . + a2nxn = b2

. . . . . . . . . . . . . . . . . . . . .an1x1 + . . . + annxn = bn

transformari elem.−→

−→

a11x1 + a11x1 + . . . + a1nxn = b1a22x2 + . . . + a2nxn = b2

. . ....

......

annxn = bn

cu aii = 0, ∀ i = 1, . . . , n .

71



2. Rezolvarea sistemului Ax = b prin metoda substitutiei ınapoi (”back-ward substitution”), adica

xn = bn

ann

xi =1

ai,i

bi − j>i

aij · x j

, j = n − 1, . . . , 1.

Exemplul 22 S˘ a se rezolve prin metoda EG sistemul Ax = b, unde

A =

2 −1 0 0−1 2 −1 00 −1 2 −10 0

−1 2

, iar b =

10100

.

Avem succesiv 2 −1 0 0 | 1

−1 2 −1 0 | 00 −1 2 −1 | 10 0 −1 2 | 0

1

2L1+L2→L2−→

2 −1 0 0 | 1

0 32 −1 0 | 1

2

0 −1 2 −1 | 10 0 −1 2 | 0

2

3L2+L3→L3−→

2 −1 0 0 | 10 3

2 −1 0 | 12

0 0 43 −1 | 4

3

0 0 −1 2 | 0

3

4L3+L4→L4−→

2 −1 0 0 | 10 3

2 −1 0 | 12

0 0 43 −1 | 4

3

0 0 0 54 | 1

.

Sistemul Ax = b de mai ınainte este 2x1 − x2 = 1

32 x2 − x3 = 1

2

43 x3 − x4 = 4

3

54 x4 = 1

si se rezolv˘ a prin substitut ie ınapoi ( x4 =4

5, x3 =

3

4(x4 +

3

4),etc).

Daca acceptam si permutari de coloane ın matricea A, un sistem de forma(7.1) se poate aduce ıntotdeauna la o forma superior triunghiulara Ax = b,dar matricea A poate fi singulara, deci etapa 2 de substitutie ınapoi nu se

mai poate aplica. Vom exemplifica acest lucru ın cele ce urmeaza.Exemplul 23 Rezolvat i prin EG sistemul

x1 + x2 + x3 + x4 = 7x1 + x2 + 2x4 = 8

2x1 + 2x2 + 3x3 = 10x1 − x2 − 2x3 + 2x4 = 0

72



1 1 1 1 | 71 1 0 2 | 8

2 2 3 0 | 10−1 −1 −2 2 | 0

L2−L1→L2L3−2L1→L3

L4+L1→L4

−→

1 1 1 1 | 7

0 0 −1 1 | 1

0 0 1 −2 | 40 0 −1 3 | 7

C 2↔C 3

−→1 1 1 1 | 7

0 -1 0 1 | 10 1 0 −2 | 40 −1 0 3 | 7

L2+L3→L3

L4−L2→L4−→

1 1 1 1 | 70 −1 1 0 | 1

0 0 -1 0 | −30 0 2 0 | 6

L4+2L3→l4−→

1 1 1 1 | 70 −1 1 0 | 10 0 −1 0 | −30 0 0 0 | 0

,

deci A = 1 1 1 1

0 −1 1 00 0 −1 00 0 0 0

; b = 7

1−30

.

Observam ca A este superior triunghiulara, dar nu si inversabila.Pentru a ınlatura ın consideratiile teoretice care urmeaza aceste situatii sia putea stabili clasa maximala de matrice pentru care putem aplica EG,introducem urmatoarea notiune.

Definitia 17 Spunem c˘ a metoda EG functioneaza pentru sistemul (7.1)dac˘ a acesta poate fi adus prin transform˘ ari elementare (7.6) la un sistem cu matricea superior triunghiular˘ a si inversabil˘ a.

Vom spunem c˘ a EG functioneaza fara pivotare dac˘ a transformarea de mai sus poate fi obt inut˘ a f˘ ar˘ a permut˘ ari de linii.

Teorema 29 EG funct ioneaza dac a si numai dac˘ a detA = 0.

Algoritm pentru EG (fara pivotare)Algoritmul pentru metoda EG este urmatorul:

1. for j = 1 to n − 1 dofor i = j + 1 to n do

m = aij/a jj ; /* a jj = 0 */aij = 0;

for k = j + 1 to n + 1 do/* coloana n + 1 contine termenul liber*/aik = aik − m ∗ a jk ;

endforendfor

endfor

73



2. xn = bn/ann;for i = n

−1 downto 1 do

sum = 0;for j = i + 1 to n do

sum = sum + aij ∗ x j ;endforxi = (1/aii) ∗ (bi − sum)

endfor

7.2 Descompunerea LU

Definitia 18 Spunem c˘ a matricea inversabil˘ a A

∈ Mn admite o

descompunere LU dac˘ a exist˘ a matricele L si U ∈ Mn, inversabile, cu L = inf.∆, si U = sup.∆ astfel ınc at A = LU.

Observatia 59 Nu orice matrice inversabil˘ a admite o descompunere LU .

De exemplu, pentru A =

0 11 1

, nu exist˘ a matricele L si U astfel ca

A = LU (scriind cele 4 ecuat ii se obt ine un sistem incompatibil).

Observatia 60 Dac˘ a exist˘ a descompunerea LU , ea nu este unic˘ a. De

exemplu, pentru A = 2 2

1 3 , dou˘ a descompuneri LU distincte sunt

L =

1 01

21

, U =

2 20 2

, respectiv

L =

2 01 2

, U =

1 10 1

.

Urmatorul rezultat stabileste o legatura importanta ıntre EG si existentaunei descompuneri LU (ce se obtine prin EG).

Teorema 30 Dac˘ a EG funct ioneaz˘ a f˘ ar˘ a pivotare, atunci matricea A ad-mite o descompunere LU .

Observatia 61 Un exercit iu simplu ne spune c˘ a matricea

L(k)−1

este

de forma

74



L(k)

−1=

1...

. ..· · · · · · 1

−mk+1,k. . .

−mk+2,k. . .

· · ·−mn,k 1

.

Exemplul 24 Determinat i o descompunere LU pentru matricea

A =

4 0 11 3 2

0 1 2

(a) folosind EG; (b) direct (prin identificare).Solut ie. (a) Ne vom folosi de demonstrat ia Teoremei 30. Pentru aceasta,efectu˘ am prima etap˘ a din algoritmul EG.

A =

4 0 1

1 3 20 1 2

L2− 1

4L1→L2−→

4 0 1

1 3 7

40 1 2

L3− 1

3L2→L3−→

4 0 1

0 37

4

0 0 1712

=

A = U ; L =

1 0 0

1

41 0

0 13 1

.

(b) A = LU =

l11 0 0l21 l22 0l31 l32 l33

· u11 u12 u13

0 u22 u23

0 0 u33

cu uii = 0, lii = 0, ∀ i = 1, 2, 3.

Sistemul obt inut este format din 9 ecuat ii si are 12 necunoscute, ceea ce ınseamn a c˘ a trebuie sa punem noi condit ii suplimentare asupra a 3 ne-cunoscute pentru a-l rezolva. Lu˘ am, de exemplu, u11 = u22 = u33 = −2.Deci avem de rezolvat sistemul

A = l11 0 0

l21 l22 0l31 l32 l33

· −2 u12 u13

0 −2 u23

0 0 −2

care ne conduce la urm˘ atoarele ecuat ii 4 = −2 l11 ⇒ l11 = −2

0 = −2 u12 ⇒ u12 = 0

75



1 = −2 u13 ⇒ u13 = −1

2

1 = −2 l21 ⇒ l21 = −12

3 = −2 l22 ⇒ l22 = −3

2

2 =1

4− 3

2u23 ⇒ u23 =

2

3(

1

4− 2) = −7

6

0 = −2 l31 ⇒ l31 = 0

1 = −2 l32 ⇒ l32 = −1

2

2 =7

12− 2 l33 ⇒ l33 =

1

2(

7

12− 2) = −17

24

Asadar, am obt inut ın acest fel

L =

−2 0 0

−1

2−3

20

0 −1

2− 7

24

, U =

−2 0 −1

2

0 −2 −7

60 0 −2

.

Observatia 62 ˆ In cazul general, pentru rezolvarea sistemului obt ¸inut la punctul (b) anterior, se pot pune condt ii asupra unora dintre necunoscute ın mai multe moduri, dup˘ a cum urmeaz˘ a:

1. Dac˘ a impunem ca lii = 1, ∀ i = 1, . . . , n, atunci descompunerea LU obt inut˘ a se numeste descompunere de tip Doolittle (asa cum se obt ine,de exemplu, prin eliminare Gauss; vezi punctul (a) de la exercit iul anterior).

2. Dac˘ a impunem uii = 1, ∀ i = 1, . . . , n, atunci descompunerea LU obt inut˘ a se numeste descompunere de tip Crout.

3. O alt˘ a variant˘ a este lii = uii, ∀ i = 1, . . . , n .

Observatia 63 Obt inerea unei descompuneri LU (atunci cˆ and este posi-bil) este foarte util˘ a dac˘ a trebuie s˘ a rezolv˘ am mai multe sisteme liniare, care

au aceeasi matrice, adic˘ a Ax = bk, k = 1, . . . , N pentru c˘ a odat˘ a obt inuta o descompunere LU pentru A, vom rezolva sistemul LUx = Ly = bk ın dou˘ a etape

I. Ly = bk care se rezolv˘ a prin substitut ie ınainte.

II. U x = y care se rezov˘ a prin substitut ie ınapoi.

76



7.2.1 Conditii suficiente pentru ca EG sa functioneze fara

pivotare

Vom prezenta ın aceasta sectiune cateva clase de matrice pentru care EGfunctioneaza fara pivotare si deci, conform Teoremei 30 exista o descom-punere LU a lui A. Mai ıntai vom introduce cateva notiuni pe care le vomfolosi ın continuare.

Definitia 19 Spunem c˘ a matricea A ∈ Mn este diagonal dominanta

dac˘ a |aii| ≥n

j=1,j=i

|aij | , ∀ i = 1, . . . , n .

Observatia 64 Dac˘ a inegalitatea de mai sus este strict˘ a, matricea se nu-meste strict diagonal dominanta.

Teorema 31 Dac˘ a matricea A este inversabil˘ a si diagonal dominant˘ a,atunci EG funct ioneaz˘ a f˘ ar˘ a pivotare.

Exemplul 25 Matricea A =

2 −1 0 0

−1 2 −1 00 −1 2 −10 0 −1 2

satisface ipotezele

Teoremei 31, deci A admite o descompunere LU .

Definitia 20 Spunem c˘ a matricea A ∈ Mn este simetrica dac˘ a A = At.

Definitia 21 Spunem c˘ a matricea A ∈ Mn este pozitiv definita daca

Ax, x > 0, ∀ x ∈ IRn \ 0, unde x, y =n

i=1xiyi.

Mentionam ın continuare doua rezultate referitoare la matrice pozitiv def-inite (vezi de ex. [?], [24]).

Propozitia 3 Matricea A ∈ Mn este pozitiv definit˘ a dac˘ a si numai dac˘ a ∆k > 0, ∀ k = 1, . . . , n, unde ∆k este determinantul submatricei obt inute prin intersect ia primelor k linii cu primele k coloane.

Propozitia 4 Dac˘ a matricea A∈ M

n este pozitiv definit˘ a, atunci aii >0, ∀ i = 1, . . . , n .

Demonstratie In Definitia 21 luand x = ei = (0,..., 0, 1, 0,...0) ∈ IRn

pentru un indice i fixat arbitrar, se obtine Aei, ei = aii > 0, ceea ceıncheie demonstratia.

77



Teorema 32 Dac˘ a matricea A ∈ Mn este simetric˘ a si pozitiv definit˘ a,atunci EG funct ioneaz˘ a f˘ ar˘ a pivotare.

Exemplul 26 Matricea A =

2 −1 0

−1. . .

. . .. . .

. . . −10 −1 2

n

satisface ipotezele Teoremei 32, deci EG funct ioneaz˘ a f˘ ar˘ a pivotare.

Urmatorul rezultat ne ofera o conditie suficienta sa existe o descompunereLU , fara a folosi metoda EG.

Teorema 33 Fie A ∈ Mn; pentru t = 1, . . . , n − 1,

not˘ am cu At matricea dat˘ a de At =

a11 a12 · · · a1t

a21 a22 · · · a2t

· · · · · · · · · · · ·at1 at2 · · · att

∈ Mt. Dac˘ a

det(At) = 0 pentru orice t = 1, . . . , n − 1, atunci exist˘ a odescompunere LU pentru A.

Demonstratie. Vom demonstra prin inductie dupa dimensiunea n ≥ 1 alui A. Pentru n = 1, afirmatia este evident adevarata (a11 = a ·b, a,b = 0).O presupunem adevarata pentru matricele de o anumita dimensiune n ≥ 1si demonstram pentru matricele de dimensiune n + 1.

Scriem pe A sub forma

A =

An cwt an+1,n+1

∈ Mn+1, (7.7)

cu c, w ∈ IRn-vectori coloana, an+1,n+1 ∈ IR si An ∈ Mn+1 ce verificaipoteza teoremei.Din ipoteza de inductie rezulta atunci ca exista Ln ∈ Mn matrice inferiortriunghiulara si inversabila si U n ∈ Mn matrice superior triunghiulara siinversabila astfel ıncat An = Ln · U n.Vrem sa gasim matricele Ln+1 ∈ Mn inferior triunghiulara si inversabila si

U n+1 ∈ Mn superior triunghiulara si inversabila cu proprietatea ca An+1 =Ln+1U n+1.

Pentru aceasta, fie Ln+1 =

Ln 0αt x

, x ∈ IR \ 0, 0, α ∈ IRn

U n+1 =

U n β 0 y

, y ∈ IR \ 0, 0, β ∈ IRn

78



Vream sa aratam ca exista x, y ∈ IR, α, β ∈ IRn astfel ıncatAn+1 = Ln+1U n+1

⇔⇔ An cwt an+1,n+1

= Ln 0αt x

· U n β 0 y

de unde obtinem Ln · U n = An si

Ln · β = cU n · α = wβ · α + x · y = an+1,n+1

Primele doua sisteme au solutie unica (Ln, U n fiind inversabile), deci α, β sunt unic determinate. Ultima ecuatie are o infinitate de solutii cu x, y = 0.Teorema este astfel complet demonstrata.

7.3 Descompunerea Choleski

Definitia 22 Spunem c˘ a A ∈ Mn admite o descompunere Choleskidac˘ a exist˘ a matricea L ∈ Mn, inf.∆ si inversabil˘ a astfel ıncˆ at A = L · Lt.

Teorema 34 Dac˘ a A ∈ Mn este simetric˘ a si pozitiv definit˘ a, atunci ea ad-mite o descompunere Choleski A = L ·Lt. Dac˘ a elementele de pe diagonala matricei L sunt stricy pozitive, atunci descompunerea este unic˘ a.

Exemplul 27 Determinat i descompunerea Choleski (dac˘ a exist˘ a) pentru

matricea A =

3 0 10 2 11 1 3

(a) folosind o descompunere LU ca ın demonstrat ia Teoremei 34;(b) direct (prin identificare)Solut ie. Matricea A este simetric˘ a si pozitiv definit˘ a (veziPropozit ia 3),deci admite o descompunere Choleski.

(a) A =

3 0 10 2 11 1 3

L3− 1

3L1→L3−→

3 0 1

0 2 10 1 8

3

L3− 1

2L2→L3−→

L3− 1

2L2→L3−→

3 0 10 2 10 0 13

16

= U , iar L =

1 0 00 1 013

12 1

Fie D1 = diagL = I 3, D2 = diag(U ) = 3 0 00 2 0

0 0 136

.

Avem A = T DS = T DT t, unde √

D :=

√

3 0 0

0√

2 0

0 0

136

;

79



Factorul Choleski pentru matricea A este

L = 1 0 0

0 1 013

12 1

·√

3 0 0

0 √ 2 00 0

136

= √

3 0 0

0 √ 2 0√ 3

3

√ 2

2

136

;

(b) Din A = L·Lt =

l11 0 0l21 l22 0l31 l32 l33

· l11 l21 l31

0 l22 l32

0 0 l33

=

3 0 10 2 11 1 3

⇒ l2

11 = 3 ⇒ l11 =√

3√ 2 · l32 = 1 ⇒ l32 = 1√

2√ 3 · l21 = 0 ⇒ l21 = 0

√ 3 · l31 = 1 ⇒ l31 = 1√

3

0 = 0l2

22 = 2 ⇒ l22 =√

2

1√ 3

· √ 3 = 1

1√ 2

· √ 2 = 1

13 + 1

2 + l233 = 3 ⇒ l2

33 = 136 ⇒ l33 =

136

si deci factorul Choleski este ın acest caz L =

√ 3 0 0

0√

2 01

√ 31

√ 2 13

6

.

7.4 Exercitii

1. Rezolvati prin EG urmatoarele sisteme

(a)

x − y + 3z = 23x − 3y + z = −1x + y = 3

(b)

x − 1

2y + z = 4

2x − y − z + t = 5x + y = 2

x − 1

2y + z + t = 5

80



(c)

x1 + x2 + x4 = 22x1 + x2 − x3 + x4 = 1

−x1 + 2x2 + 3x3 − x4 = 43x1 − x2 − x3 + 2x4 = −3

2. Rezolvati sistemul de la exercitiul 1, determinand mai ıntai o descom-punere LU pentru matricea A.

3. Determinati care din urmatoarele matrice sunt (i) simetrice ; (ii)pozitiv definite; (iii) strict diagonal dominante; (iv) inversabile.

(a)

2 11 3

; (b)

-2 11 -3

; (c)

2 1 00 3 21 2 4

.

4. Pentru A = α 1 0−β 2 10 1 2

, gasiti α si β reale pentru care

(a) A este strict diagonal dominanta.

(b) A este simetrica.

(c) A este pozitiv definita.

5. Pentru matricea de la exercitiul 4, pentru α = 2 si β = 1 determinatio descompunere Choleski.

81



Capitolul 8

Aproximarea valorilor si

vectorilor proprii

Definitia 23 Pentru o matrice A ∈ Mn, num˘ arul λ ∈ C se numeste valoare proprie dac˘ a este r˘ ad˘ acin˘ a a polinomului caracteristic, adic˘ a

P A(λ) = det(A − λI n) = 0, (8.1)

iar x ∈ Cn \ 0 se numeste vector propriu corespunz˘ ator valorii proprii λ dac˘ a

Ax = λx. (8.2)

Vom nota cu σ (A) spectrul matricei A∈ Mn

, adica multimea valorilorproprii, iar prin ρ (A) raza spectrala a acesteia definita prin

ρ(A) = max|λ|, λ ∈ σ(A). (8.3)

Definitia 24 O matrice A ∈ Mn se numeste diagonalizabila dac˘ a ea este asemenea cu o matrice diagonal˘ a D ∈ Mn, adic˘ a exist˘ a P ∈ Mn

inversabil˘ a astfel ıncˆ at

P −1AP = D = diag(λ1,...,λn), (8.4)

si respectiv triangularizabila dac˘ a este asemenea cu o matrice U ∈ Mn

superior triunghiular˘ a, adic˘ a exist˘ a S ∈ Mn inversabil˘ a astfel ıncˆ at

S −1AS = U. (8.5)

Vom nota aceste relatii prin A ≈ D, respectiv A ≈ U .

82



de rotatie Givens definita prin (vezi de exemplu [23])

Ω =

p q 1

. . .

1cos θ sin θ

1. . .

− sin θ cos θ1

. . .

1

p

q

(8.9)

Deoarece Q este ortogonala, daca definim matricea B ∈ Mn prin

B = ΩtAΩ (8.10)

atunci avem (vezi Exercitiul 17 de la Cap. 8)

B2F =

ni,j=1

b2ij =

ni,j=1

a2ij = A2

F . (8.11)

Lema 4 (i) Dac˘ a a pq = 0, atunci exist˘ a si este unic θ ∈ π4 , π4 \ 0 astfel ınc at

bqp = 0. (8.12)

(ii) Pentru θ de la (i) avem

ni=1

b2ii =

ni=1

a2ii + 2a2

pq. (8.13)

Demonstratie. Evidenta. Metoda Jacobi pentru aproximarea valorilorsi vectorilor proprii ai matricei simetrice A este urmatoarea: se definesteA1 = A si se genereaza sirul de matrice (Ak)k≥1 prin

Ak+1 = ΩtkAkΩk, k ≥ 1, (8.14)

unde Ωk este de forma (8.9) cu θ = θk dat de (vezi (??))

ctg2θk =ak

qq − ak pp

ak pq

, (8.15)

84



(prin akij am notat elementele matricei Ak, k ≥ 1). Indicii p si q se aleg (ın

cazul metodei Jacobi clasice) prin

|ak pq| = max

i= j|ak

i,j |. (8.16)

Observatia 65 La fiecare iterat ie k ≥ 0, ın matricea Ak corespunz˘ atoare anul˘ am elementul de pe pozit ia ( p, q ), alta decˆ at cea de la pasul anterior.Totusi, pe parcursul calculelor s-ar putea ca pe unele pozit ¸ii (deja anulate)s˘ a apar˘ a din nou elemente nenule, acestea ıns˘ a au va-lori (absolute) din ce ın ce mai mici.

In ipotezele si cu constructiile anterioare sunt adevarate urmatoarele rezul-

tate (pentru demonstratii vezi [23]).Teorema 35 (Convergenta valorilor proprii)Sirul (Ak)k≥1 dat de metoda Jacobi clasic˘ a (8.14) − (8.16) este convergent la matricea

D = diag(λπ(1),...,λπ(n)) (8.17)

unde π este o permutare a mult imii 1,...,n, iar σ(A) = λ1,...,λn este spectrul lui A.

Fie acum Qk ∈ Mn matricea ortogonala

Qk = Ω1Ω2...Ωk, k ≥ 1 (8.18)

cu Ωk din (8.14)-(8.15).

Teorema 36 (Convergenta vectorilor proprii)Dac˘ a spectrul lui A cont ine valori proprii mutual distincte, atunci sirul (Qk)k≥1 din (8.18) converge la o matrice ortogonal˘ a ale c˘ arei coloane formeaz˘ a un sistem ortonormat de vectori proprii ai lui A.

Observatia 66 Metoda Jacobi se refer˘ a doar la aproximarea valorilor si vectorilor proprii pentru o matrice simetric˘ a. Pentru matrice mai generale,se pot consulta lucr˘ arile [ 3 ], [ 5 ] , [ 6 ] si [ 17 ] .

8.2 Metoda puterilor

Fie A ∈ Mn matrice cu proprietatea ca exista o singura valoare proprie demodul maxim si exista o baza de vectori proprii, adica σ(A) = λ1,...,λn,

|λ1| > |λ2| ≥ ... ≥ |λn| (8.19)

85



si u1,...,un o baza ın Cn astfel ıncat

Aui = λiui, i = 1, . . . , n (8.20)

Pentru k ≥ 1 arbitrar fixat si x0 ∈ Cn de forma,

x0 = a1u1 + a2u2 + ... + anun, (8.21)

cu ai ∈ IR, i = 1, . . . , n, a1 = 0, definim

x1 = Ax0, x2 = Ax1, . . . , xk = Axk−1, (8.22)

decixk = Akx0. (8.23)

Deoarece a1

= 0, putem presupune fara a restrange generalitatea ca

x0 = u1 + a2u2 + ... + anun. (8.24)

Din (8.23) si (8.20) rezulta egalitatea

xk = λk1u1 + λk

2a2u2 + · · · + λknanun, (8.25)

care poate fi rescrisa sub forma

xk = λk1

u1 +

λ2

λ1

k

a2u2 + · · · +

λn

λ1

k

anun

(8.26)

Cum λ j

λ1 < 1, ∀ j = 2, . . . , n (vezi (8.19)), rezulta λ j

λ1k

k→∞−→ 0, pentru

orice j = 2, . . . , n. Din (8.26) obtinem ca xk este de forma

xk = λk1[u1 + εk], (8.27)

unde εk → 0, k → ∞. Fie acum Φ o functionala liniara si continua pe Cn

astfel ıncat Φ(u1) = 0. Din (8.27) rezulta Φ(xn) = λk1[Φ(u1) + Φ(εn)], deci

rk :=Φ(xk+1)

Φ(xk)= λ1

Φ(u1) + Φ(εk+1)

Φ(u1) + Φ(εk)

→ λ1

adica λ1 este aproximat prin elementele sirului (rk)k

≥1 si demonstratia este

completa.Algoritmul de calcul pentru metoda puterilor este urmatorul.

Metoda puterilor

86



Introduceti n,A,x, N si Φ o functionala linara si continuafor k = 1, . . . , N do

y = Axr = Φ(y)/Φ(x)x = y/yAfiseaza k,x,r

endforAfiseaza x.

Observatia 67 Ultima valoare afisat˘ a pentru r va fi aproximat ia pentru λ1, iar ultima valoare pentru x va fi aproximat ia pentru vectorul propriu corespunz˘ ator valorii proprii λ1.

Exemplul 28 Fie A = 6 5−

52 6 −22 5 −1

, x = (−1, 1, 1), N = 28 si Φ = x2.

Algoritmul dat de metoda puterilor va afisa urm˘ atoarele valori:

k = 0 x0 = (−1.00000, 1.000000, 1.0000000)k = 1 x1 = (−1.00000, 0.333333, 0.333333) r0 = 2.0k = 2 x2 = (−1.00000, −0.111111, −0.111111) r1 = −2.0k = 3 x3 = (−1.00000, −0.407407, −0.407407) r2 = 22.0. . . . . . . . .k = 6 x6 = (−1.00000, −0.824417, −0.824417) r5 = 6.71508. . . . . . . . .k = 28 x28 = (

−1.00000,

−0.999977,

−0.999977) r27 = 6.00007

Pentru matricea A de mai ınainte, raza spectral˘ a este 6, iar un vector propriu asociat (1, 1, 1)t.

8.3 Metoda puterilor inverse

Pentru calculul celei mai mici valori proprii ın modul, se foloseste metodaputerilor inverse. Astfel, fie A inversabila si σ(A) = λ1,...,λn cu

|λ1| ≥ |λ2| ≥ ... ≥ |λn−1| > |λn| > 0. (8.28)

|λ−1n | > |λ−1

n−1| ≥ ... ≥ |λ−11 | > 0 (8.29)

si se foloseste ın acest caz metoda puterilor pentru calculul lui λ−1n , aplicata

matricei A−1, pentru care σ(A−1) =

1

λ1,...,

1

λn

.

Nu se recomanda calculul lui A−1, ci ın loc de xk+1 = A−1xk, se rezolva

87



sistemul Axk+1 = xk. Aceasta poate fi facuta ın mod eficient folosindmetoda bazata pe descompunerea LU a lui A.

8.4 Exercitii

1. Sa se arate ca doua matrice asemenea au acelasi spectru.

2. (Teorema lui Gershgorin) Pentru A ∈ Mn, σ(A) ⊆n

i=1Di, unde

Di =

z ∈ C/ |z − aii| ≤

n j=1,j=i

|aij|

.

(Aceasta teorema se mai numeste si teorema de localizare a valorilorproprii)

3. Sa se arate ca orice matrice strict diagonal dominanta este inversabila.

4. Fie A matrice simetrica si pozitiv definita. Atunci I −A este simetricasi pozitiv definita ⇔ ρ(A) < 1.

5. Pentru A =

2 0 −1−2 −10 0−1 −1 4

determinati o valoare aproximativa

pentru ρ(A) efectuand doi pasi ai metodei puterilor cu x0 = (1, 1, 1)t

si Φ(x) = x∞.

6. Determinati complexitatea aritmetica a algoritmului de la metodaputerilor.

88



Bibliografie

[1] Bakhvalov N., Methodes numeriques , Editions MIR, Moscou 1976.

[2] Branzanescu V., Stanasila O., Matematici speciale - partea I , Ti-pografia Institului Politehnic Bucuresti, 1985.

[3] Bourden R. L., Faires J. D., Reynolds A. C., Numerical Analysis -second edition , Prindle, Weber and Schmidt, Boston, Massachusetts,1981.

[4] Ciarlet P.G., The finite element method for elliptic problems , North-Holland, 1979.

[5] Ciarlet P.G., Introduction a l’analyse numerique matricielle et a l’optimisation , Masson, Paris,1982.

[6] Demidovici B.P., Maron I.A., Computational Mathematics , MIR Pub-lishers, Moscow, 1981.

[7] Dragos L., Metode matematice ın aerodinamic˘ a , Editura AcademieiRomane, Bucuresti, 2000.

[8] Dragos L., Popescu M, Certain quadrature formulae of interest in aero-dynamics , Rev. Roum. Math. Pures Appl., XXXVII(7)(1992), pp.587-599.

[9] Gantmacher F.R., The theory of matrices (vol. I si II), Chelsea Publ.Comp., New York 1959.

[10] Golub G. H., van Loan C. F., Matrix computations - third edition , The

John’s Hopkins Univ. Press, Baltimore,1996.

[11] Hackbusch W., Elliptic differential equations - Theory and numerical treatment , Springer-Verlag, Berlin Heidelberg 1992.

[12] Henrici P., Elements of numerical analysis , John Willey and Sons,Nwy York 1964.

89



[13] Householder A. S., The theory of matrices in numerical analysis , Blais-dell Publ. Comp., New York, 1964.

[14] Juncu Gh., Popa C., Introducere ın metoda multigrid - aplicat ii pe calculator , Editura Tehnica, Bucuresti 1991.

[15] Kantorovici L.V., Akilov G.P., Analiz˘ a funct ional˘ a (trad. din limba rus˘ a), Ed. St. si Encicl., Bucuresti 1986.

[16] Kinkaid D.R., Cheney W., Numerical Analysis: Mathematics of Sci-entific Computing , Brooks/Cole Publishing Company, Pacific Groove,California 1991.

[17] Meyer C.D., Matrix analysis and applied linear algebra , SIAMPhiladelphia 2000.

[18] Micula Gh., Funct ii spline si aplicat ii , Editura Tehnica, Bucuresti1978.

[19] Overton M.L., Numerical computing with IEEE floating point arith-metic , SIAM Philadelphia 2001.

[20] Pelican E., Analiz˘ a numeric˘ a - complemente, exercit ii, programe de calcul , va apare ın Editura MatrixRom, Bucuresti, 2006.

[21] Popa C., Pelican E., Introducere ın analiza numeric˘ a , Editura Ma-trixRom, Bucuresti, 2005.

[22] Popa C., Introducere ın analiza numeric˘ a , Editura EUROBIT,Timisoara 1996.

[23] Popa C., Analiz˘ a numeric˘ a matriceal˘ a , Editura EUROBIT, Timisoara1996.

[24] Popa C. si colectiv., Analiza numeric˘ a. Complemente. Exercit ii. Pro-grame de calcul , Tipografia Universitatii Ovidius, Constanta 1996.

[25] Popoviciu T., Analiza numeric˘ a - not ¸iuni introductive de calcul aprox-imativ , Ed. Academiei Romane, Bucuresti 1991.

[26] Siretchi Gh., Calcul diferent ial si integral, vol. I si II , Ed. St. si Encicl.,Bucuresti 1985.

[27] Suli E., Mayers D., An introduction to numerical analysis , CambridgeUniversity Press, 2003.

[28] Varga R., Matrix iterative analysis , Prentice Hall, New York,1962.

90



[29] Weiss R., Parameter-free iterative linear solvers , Mathematical Re-search Series, 97, Akademie Verlag, Berlin, 1996.

[30] Young D. M., Iterative solution of large linear systems , AcademicPress, New York,1971.

Analiza numerica.pdf curs de matematica

Documents

Transcript of Analiza numerica.pdf curs de matematica