PARTEA I

PRINCIPIILE MATEMATICII NUMERICE

ANALIZA STABILITĂŢII ŞI ERORILOR

1.1. Probleme stabile (sau bine puse)

Considerăm problema: calculează x astfel încât

0dxF , (1)

unde d este mulţimea datelor de care depinde soluţia x şi F este relaţia funcţională între x şi d.

În funcţie de problema particulară rezolvată, x şi d pot fi din R, din Rn (vectori) sau din R

mxn

(matrice).

Definiţia 1. Dacă F şi d sunt date (cunoscute), problema (1) se numeşte directă. Dacă F

şi x sunt date, problema (1) se numeşte inversă. Dacă d şi x sunt date, problema (1) se numeşte

problemă de identificare.

Definiţia 2. Problem (1) se numeşte stabilă sau bine pusă dacă soluţia x este unică şi

depinde continuu de datele d. Dacă problema (1) nu este stabilă, pentru rezolvarea sa numerică

este necesară transformarea ei într-o problemă stabilă.

Observaţie. Prin dependenţă continuă de date înţelegem faptul că perturbaţii mici ale

datelor d determină modificări “mici” ale soluţiei x. Cu alte cuvinte, dacă d este perturbarea

adusă datelor şi x este perturbarea indusă soluţiei şi are loc relaţia

0 ddxxF , (2)

atunci,

dK , 0 astfel încât dxd

ddKxd ,

unde x

. este norma pe spaţiul soluţiilor şi d

. este norma pe spaţiul datelor.

În continuare, pentru simplificare, vom nota ambele norme cu . , faptul că o normă

este definită pe un spaţiu sau celălalt fiind înţeles din context.

Exemplul 1. Fie polinomul 112 24 aaxaxxp . Problema 0, axF ,

unde xpaxF , nu este stabilă deoarece există o variaţie discontinuă a numărului

rădăcinilor din R ale lui xp funcţie de parametrul continuu a.

Evident, notând 02 xt , rezultă

01122 aatat ,

cu soluțiile 11 at și at 2

Dacă 1a sunt 4 soluţii distincte în R

Dacă 1a , 1,0 21 tt şi sunt 4 soluţii în R (două confundate)

Dacă 10 a sunt 2 soluţii distincte în R

Dacă 0a sunt 2 soluţii confundate în R

Altfel, problema nu are soluţii în R

Observaţie. Condiţia de unicitate a soluţiei poate fi îndeplinită, dacă prin soluţie

înţelegem un vector cu fiecare componentă rădăcină a polinomului (vezi exemplul 2). Pentru

ca acest lucru să fie posibil trebuie ca, indiferent de date (aşa cum este formulată problema,

a R), numărul rădăcinilor să rămână acelaşi. De exemplu, problema devine este stabilă

pentru mulţimea datelor satisfăcând 4a .

Definiţia 3. Pentru problema (1) indicatorul de condiţie relativ este definit prin

dd

xxdK

Dd

sup

unde D este o vecinătate a originii şi desemnează setul perturbaţiilor admisibile asupra

datelor d pentru care problema perturbată (2) are sens. Dacă 0d sau 0x este definit

indicatorul de condiţie absolut, prin

d

xdK

Ddabs

sup

Definiţia 4. Problema (1) este bine condiţionată dacă dK este “mic” pentru orice

perturbaţie admisibilă.

Observaţie. Chiar dacă indicatorul de condiţie formal este infinit, problema poate fi

reformulată echivalent astfel încât să devină problemă stabilă (vezi exemplu 2).

Dacă problema (1) admite o soluţie unică, atunci există o funcţie G, numită funcţie care

rezolvă problema (“rezolvitor”), definită pe mulţimea valorilor datelor d cu valori în mulţimea

valorilor soluţiei şi astfel încât

0,, ddGFdGx (3)

Pe baza relaţiei(2), rezultă

ddGxx (4)

Dacă G este diferenţiabilă, din dezvoltarea în serie Taylor, prin aproximare de ordinul I,

obţinem,

dOddGdGddG ' , pentru 0d

ddGdGddG '

Pe baza relaţiilor (3) şi (4), obţinem,

dGddGxddGx

dd

dGddG

dd

dGdGddG

dd

xxdK

DdDdDd

'supsupsup

Observaţie. Dacă :G Rn R

m, atunci dG' este matricea Jacobian calculată în

vectorul d.

Deoarece, dacă . este norma unui vector, rezultă că funcţia x

AxA

x 0

sup

este norma

indusă pe mulţimea matricelor, obţinem,

dG

d

ddG

Dd

''

sup

şi deci

dG

ddGdK '

Similar,

dGdKabs '

Observaţie. Dacă :G Rn R

m, atunci dG' este matricea Jacobian calculată în

vectorul d.

Exemplul 2. Fie ecuaţia algebrică de gradul II,

0122 pxx , cu 1p

Deoarece 14 2 p , rădăcinile sunt 12 ppx .

În acest caz 12, 2 pxxpxF , pd şi x este vectorul Txx , . Funcţia

“rezolvitor” este,

:G RR2, TxxpG ,

TpppppG 1,1 22

Obţinem,

11

11

'

2

2

p

p

p

p

pG

Considerând 22

2,.. yx

y

x

, rezultă,

2411 22

22

2 ppppppG

1

24

11

11'

2

22

2

2

2

p

p

p

p

p

ppG

pp

Obţinem,

1

'2

p

p

pG

ppGpK

Concluzii

1. În cazul rădăcinilor bine separate (de exemplu 2p ), problema 0, pxF este

bine condiţionată.

2. În cazul rădăcinilor multiple, pentru 1p , funcţia :G RR2,

TpppppG 1,1 22 nu este diferenţiabilă, deci pK nu poate fi calculat în

funcţie de G.

În plus, dacă 1p , dar este apropiat de 1, pK are valori mari şi problema nu este

bine condiţionată. Problema poate fi reformulată într-o manieră echivalentă şi astfel încât să

devină stabilă.

Fie 11

,2

2

xt

txtxF , 0, txF , cu 12 ppt .

Evident, pt

t2

1 2

.

Cele două rădăcini sunt t

x1

, tx

1

11

2

2

pppp şi coincid pentru

1 pt . Pentru funcţia considerată, rezultă

T

tttG

1

T

ttG

2

11'

t

t

tttG

11 4

2

2

2

4

4

111'

t

t

ttG

Rezultă,

1' tG

ttGtK

deci problema este stabilă în această reprezentare.

1.2. Stabilitatea şi convergenţa metodelor numerice

Considerăm că problema (1) este stabilă. O metodă numerică pentru aproximarea

soluţiei problemei (1) este în general o secvenţă de probleme ce depind de 1n ,

1,0, ndxF nnn (5)

şi astfel încât xxn când n , adică soluţia numerică este convergentă către soluţia

exactă. Pentru aceasta trebuie ca ddn şi nF să “aproximeze” F când n .

Definiţia 5. Dacă setul de date d din (1) este admisibil pentru nF , problema numerică

(5) se numeşte consistentă dacă

0,,, dxFdxFdxF nn când n . (x este soluţia problemei 1)

O metodă se numeşte consistentă în sens tare dacă

0,,, dxFdxFdxF nn pentru orice n (nu numai dacă n ).

Definiţia 6. În cazul metodelor iterative, problema (5) are forma,

qndxxxF nqnnnn ,0,,...,, 1 (6)

cu 110 ,...,, qxxx valori date. În acest caz consistenţa în sens tare revine la verificarea

proprietăţii,

qndxxxFn ,0,,...,, (7)

Exemplul 4. Fie metoda iterativă de aproximare a rădăcinii unice a funcţiei

:f RR, dată de relaţia

1,' 1

11

n

xf

xfxx

n

nnn , 0x dat. (metoda Newton)

Conform definiţiei, metoda este consistentă în sens tare, deoarece

1,0,, nfFn , unde

1,'

,,1

111

n

xf

xfxxfxxF

n

nnnnnn

Într-adevăr,

1,0'

,, nf

ffFn

, deoarece este rădăcina funcţiei f.

Definiţia 7. O metodă numerică asociată problemei (1) este stabilă (bine pusă) dacă,

pentru orice 1n să aibă loc proprietatea,

nnnnn ddKxddK ,:,0

Similar problemei (1), pentru fiecare problemă din secvenţa (5) sunt definiţi indicatorii

nn

nn

Ddnn

dd

xxdK

nn

sup

n

n

Ddnnabs

d

xdK

nn

sup,

şi definim:

indicatorul condiţie asimptotică relativ: nnknk

nnum dKdK

suplim

indicatorul condiţie asimptotică absolut: nnabsknk

nnumabs dKdK ,suplim

Metoda numerică este bine condiţionată dacă numK este suficient de “mic” pentru toate

datele admisibile nd .

Similar problemei (1), dacă presupunem existenţa câte unei soluţii unice nx a celei de-a

n-a probleme din secvenţa (5), atunci există câte o funcţie care rezolvă problema (“rezolvitor”),

nG definită prin,

1,0,, nddGFdGx nnnnnnn

Dacă nG este diferenţiabilă, similar indicilor condiţie definiţi pentru problema (1), sunt

calculaţi indicii condiţie asimptotică prin,

nn

nnnnn

dG

ddGdK '

nnnnabs dGdK ',

Exemplul 5. Fie problema din exemplul 2 şi considerăm 1p (cazul rădăcinilor

separate). Dacă 12 ppx este evaluată direct algoritmul nu este stabil. Într-adevăr,

pentru p cu valori mari, calculul formulei

12 pp

0 deapropiat este12pp induce erori mari de calcul datorită

reprezentării aritmetice finite (vezi §1.5, situaţia subdepăşirii de domeniu). Variantele de

rezolvare a problemei de subdepăşire:

calculează 12 ppx şi, evident,

x

x1

rezolvă ecuaţia pentru domeniu ce conţine numai valoare x cu metoda Newton

1,' 1

1

1

nxf

xfxx

n

n

nn , 0x dat

1,22

12

1

1

2

11

npG

px

pxxxx n

n

nnnn

px

pxxxxpxxF

n

nnnnnn

22

12,,

1

1

2

111

px

xpG

n

nn

1

2

1

2

1

2

1

2

1

2

1'

px

xpG

n

n

n

px

p

pG

ppGpK

nn

nn

1

'

px

ppK

nknk

num

1

suplim

Dar 12

1 ppx

nn deoarece 12 ppxxn sau

12 ppxxn

Rezultă,

12

p

ppKpK num

n

deci metoda nu este bine condiţionată când 1p , altfel metoda este bine condiţionată.

Definiţia 8. Metoda numerică (5) este convergentă dacă,

nnn ddxdxndnn

nn

,,

:,,,0

00

00

unde d este setul datelor admisibile pentru (1), dx este soluţia problemei (1) şi nn ddx

este soluţia problemei (5) cu setul de date ndd .

Observaţie. Pentru a verifica proprietatea de convergenţă a metodei numerice (5) este

suficientă verificarea următoarei proprietăţi,

2

nnn ddxddx

Într-adevăr, din regula triunghiului obţinem,

2 def.2

,0

nnnn

nnnn

nnnnnn

dnKddxddxx

ddxddxddxdx

ddxddxddxdxddxdx

Prin alegerea lui nd astfel încât 2

,0

ndnK , rezultă,

nn ddxdx

Cele mai utilizate măsuri ale convergenţei şirului 1nnx sunt

eroarea absolută: nn xxxE

eroarea relativă: 0

xx

xxxE

n

nrel

dacă soluţia este vector sau matrice sunt utilizate erorile pe componente,

ji

jin

jin

c

relx

xxxE

,

,

,max

Observaţii.

1. Dacă problema (1) este bine pusă, pentru ca metoda numerică să conveargă ea

trebuie să fie stabilă.

2. Pentru orice metodă numerică cu proprietatea că este consistentă, stabilitatea este

echivalentă cu convergenţa (teorema Lax-Richtmyer)

1.3. Analiza a priori şi a posteriori a stabilităţii

Analiza stabilităţii unei metode numerice poate fi realizată a priori în una din

următoarele variante

analiza înainte (forward analysis), care furnizează o limită a variaţiei soluţiei,

nx , datorate perturbaţiei în date şi erorile intrinseci ale metodei numerice;

analiza înapoi (backward analysis), care are ca scop estimarea perturbaţiilor care

sunt aduse datelor unei probleme specificate pentru a obţine rezultatul calculat

de metodă. Cu alte cuvinte, dată fiind soluţia calculată nx̂ , scopul este de a

determina perturbaţia în setul de date, nd , astfel încât 0,ˆ nnn ddxF .

Analiza a priori poate fi utilizată şi pentru investigarea proprietăţii de convergenţă a

metodelor numerice, caz în care este referită drept analiza a priori a erorilor.

Analiza a posteriori are ca scop evaluarea erorii nxx ˆ ca funcţie a reziduurilor,

dxFr nn ,ˆ

prin intermediul unor constante numite factori de stabilitate, unde nx̂ este soluţia calculată

numeric, aproximare a lui x, soluţia problemei (1).

1.4. Surse de eroare în modelele de calcul numeric

Dacă problema numerică (5) este o aproximare a problemei matematice (1) care, la

rândul ei este modelul unei probleme reale (fizice) PP, problema (5) este numită model de

calcul al problemei PP.

Eroarea globală, e, este definită ca diferenţa dintre soluţia calculată de modelul de

clacul (5), nx̂ , şi soluţia problemei PP, PPx modelate prin problema (1) care are soluţia x.

Eroarea globală poate fi interpretată ca cm eee , unde,

me , eroarea modelului matematic, PPm xxe

ce , eroarea modelului de calcul, xxe nc ˆ

Similar, eroarea modelului de calcul este exprimată prin suma următoarelor valori,

ae , eroarea indusă de algoritmul numeric şi erorile de calcul datorate

reprezentărilor numerelor reale în calculator (ca mulţime finită de numere cu un

număr finit de zecimale)

ne , eroarea determinată de procesul de transformare a problemei într-una

discretă, xxe nn

Situaţiile desrise mai sus sunt sintetizate în figura 1.1.

PP: PPx

nx̂

me e

ce

ne ae

Figura 1.1

0, dxF

0, nnn dxF

În general, sursele erorilor sunt următoarele,

erori datorate modelului, care pot fi corectate printr-o alegere potrivită a

modelului matematic;

erori în date, care pot fi controlate prin creşterea acurateței măsurătorilor;

erori de trunchiere, care provin din limitarea numărului de iteraţii (etape) ale

algoritmului de calcul al unei soluţii;

erori de rotunjire, care sunt cauzate de faptul că în calculator pot fi reprezentate

doar un număr finit de numere reale şi acestea au un număr finit de cifre

zecimale.

1.5. Reprezentarea numerelor în calculator. Erori de rotunjire

Orice operaţie pe calculator este afectată de erori de rotunjire (rounding errors sau

roundoff), datorate faptului că din mulţimea numerelor reale R, doar o mulţime finită pot fi

reprezentate în calculator.

Reprezentarea numerelor reale (sistemul poziţional)

Fie o bază, , , x cu un număr finit de cifre xk, , k=-m, …, n

( ). Notaţia convenţională:

(8)

este numită reprezentarea poziţională a lui x în raport cu baza β , unde s este semnul iar

punctul dintre x0 şi x-1 este punctul zecimal (dacă β este 10) sau punctul binar (dacă β este 2).

Reprezentarea (8) este echivalentă cu

.

Exemplul 6.

x10 = 425.31

x6 = 425.31

Un număr raţional poate avea un număr infinit de cifre în reprezentarea într-o bază şi un

număr finit de cifre în reprezentarea într-o altă bază. De exemplu 1/3 în baza 10 este

reprezentat ca 0.(3) iar în baza 3 reprezentarea este 0.1.

Orice număr real poate fi aproximat prin numere cu reprezentare finită. Pentru

bază finită, este îndeplinită următoarea proprietate:

(9)

unde are un număr finit sau infinit de cifre.

Fie (cu număr finit sau infinit de cifre) şi r ,

oarecare. Atunci numerele

şi

au cîte r cifre şi sînt astfel încît

şi =

Dacă r este astfel încît , atunci rezultă (9) cu sau .

Reprezentarea numerelor reale în calculator în sistem flotant (virgulă mobilă)

Dacă sînt N poziţii de memorie (biţi) alocaţi reprezentării unui număr n, atunci ei sînt

repartizaţi astfel: unul pentru semn, N-k-1 pentru partea întreagă şi k pentru partea fracţionară.

, (10)

(sistem virgulă fixă). (10) este echivalent cu

(11)

Sistemul virgulă fixă limitează valorile minime şi maxime care pot fi reprezentate, cu

excepţia situaţiilor în care N este foarte mare. Dezavantajul poate fi depăşit prin utilizarea unei

scale variabile pentru (11), obținând astfel reprezentarea în sistemul flotant sau virgulă mobilă,

astfel:

(în general şi 1 tm

obținându-se forma normalizată),

unde:

este numărul cifrelor semnificative ,

este mantisa,

este exponentul.

Fie e: (în general şi ). Cele N poziţii disponibile sînt distribuite

astfel:

1 pentru semn

t pentru mantisă

N-t-1 pentru exponent.

Pentru reprezentarea în simplă precizie se folosesc 32 biţi iar pentru reprezentarea în

dublă precizie se folosesc 64 biţi.

Fig. 1.2

Fie , . Evident, dacă

=> .

În plus, , , .

.

Cel mai mic număr pozitiv din normalizat este . Cel

mai mare număr pozitiv din este (

.

Observaţie. Pentru ,2

Dacă F nu e normalizat, atunci .

Exemplul 7. Setul conţine, în reprezentarea normalizată, următoarele

numere pozitive:

e=2:

e=1:

e=0:

e=-1:

exponent mantisă

exponent mantisă

1 8 23

1 11 52

semn

semn

În varianta denormalizată sînt adăugate numere pentru e=-1 (de exemplu, pentru e=2,

, deci există deja reprezentarea):

Observaţii: 1. În mulţimea numerele nu sînt egal distanţate, ele devenind mai dense către cel mai

mic număr care poate fi reprezentat.

2. Dacă , consecutive, cea mai mică distanţă posibilă este

iar cea mai mare este , unde este denumită epsilon calculată şi desemnează

distanţa între 1 şi cel mai apropiat minor din F, deci este cel mai mic număr din F cu

proprietatea că .

3. Pe intervale de tipul numerele din F sînt egal distanţate de .

4. Spre deosebire de distanţa absolută dintre două numere consecutive din F (dea descrisă la 2. şi

3.), distanţa relativă are un comportament periodic şi depinde de mantisa m.

Fie . Numărul următor se află la distanţa , deci

distanţa relativă este .

În intervalul , descreşte cu creşterea lui x (în reprezentarea normalizată

mantisa variază între şi exclusiv). Cînd x devine distanţa relativă revine la

valoarea şi reîncepe descreşterea. Fenomenul oscilatoriu este cu atît mai pronunţat cu cît

β este mai mare. Este un motiv în plus pentru care este preferat β mic.

Standardul folosit este IEC559:

simplă precizie: ,

dublă precizie:

cu includerea reprezentării denormalizate

format extins.

Rotunjirea numerelor reale în operaţii pe F

Fie (parametrii sînt fixaţi) reprezentarea numerelor reale în calculator.

Prima problemă de ordin practic este reprezentarea în F a oricărui număr real (chiar dacă

, nu rezultă că , op fiind un operator binar). În consecinţă, este necesară

definirea unei aritmetici pe F. Cea mai simplă abordare este rotunjirea numerelor. Fie în

reprezentarea normalizată (12); x este substituit cu reprezentarea ,

, unde .

Evident, 1) dacă =>

2) dacă , cu => .

Observaţii:

1. dacă , fl(x) nu este definit; dacă x este rezultatul unei operaţii se

află în situaţia depăşirii de domeniu (overflow).

2. Dacă , operaţia de rotunjire poate fi definită; dacă x este rezultatul unei

operaţii pe F se află în situaţia de subdepăşire de domeniu (underflow).

3. Dacă eroarea absolută şi eroarea relativă rezultate prin

substituirea lui x cu fl(x) este stabilită prin următorul rezultat:

, cu ,

unde ( este precizia calculatorului sau unitatea roundoff).

Obţinem:

pentru eroarea relativă ,

pentru eroarea absolută

tttt

te aaaaaaaaxflxxE ~............. 212121

. Dar

1

2121212

~.............

tttt aaaaaaaa ,

deci

.