Met Ode Nu Me Rice

8/14/2019 Met Ode Nu Me Rice

1/313

FitVisible
http://view/


2/313

Aceasta este versiunea electronic a a c art ii Metode Numerice publicata deEditura Tehnic a. Cartea a fost culeasa folosind sistemul LATEX a lui LeslieLamport, o extindere a programului T EX scris de Donald Knuth. Versiuneaelectronic a foloseste formatul Portable Document Format (PDF) elaborat deAdobe Systems . Traducerea formatului L ATEX n PDF a fost realizat a cu pro-gramul pdflatex elaborat de Han The Thanh. Hiperleg aturile din versiuneaelectronic a au fost generate automat folosind pachetul hyperref al lui SebastianRahtz.

COPYRIGHT c 1998, Corneliu Berbente, Sorin Mitran, SilviuZancu

Toate drepturile asupra edit iei electronice sunt rezervate autorilor. Nu estepermis a tip arirea cont inutului acestei edit ii far a consimt amantul scris al auto-rilor.

COPYRIGHT c 1997, Editura TehnicaToate drepturile asupra edit iei tip arite sunt rezervate editurii.

Adresa: EDITURA TEHNIC APiat a Presei Libere, 133 Bucuresti, Rom aniacod 71341

Redactor: ing. Vasile BuzatuTehnoredactor: Diana JilavuCoperta: Sorin Mitran

Bun de tipar: 15.11.1997; Coli tipo: 17,75CZU: 519.6ISBN 973-31-1135-X
http://www.adobe.com/


3/313

PREFAT A

Extraordinara dezvoltare a tehnicii de calcul din ultima perioad a permite si-mularea unui tot mai mare num ar de probleme zice, ingineresti sau economice.In paralel, a avut loc o dezvoltare a programelor disponibile zicianului, ingine-rului sau economistului, oferindu-le o bogat a gam a de algoritmi pentru rezolvarea unor aplicat ii concrete. Utilizarea acestei bog at ii de tehnici si informat ii necesit a ns a, o baz a teoretic a solid a pentru a efectiv folosite.

Reprezentarea printr-un num ar nit de cifre semnicative a numerelor n calculator introduce dicult at i extrem de mari n asigurarea condit iilor pentru aplicarea unora din not iunile fundamentale ale matematicilor moderne, legatede procesul de trecere la limit a, amend and astfel utilizarea ecient a a unor te-oreme de baz a din analiz a. In schimb, se introduc erorile de rotunjire a c aror propagare, n interact ie cu alte tipuri de erori (inerente sau de metod a) estegreu de urm arit. Prinre consecint e, se poate nt ampla ca varainate echivalenteteoretic (spre exemplu pe baza unor teoreme privind unicitatea solut iei) s a duc a,numeric, la rezultate foarte diferite. Ca urmare, este explicabil a tendint a de a sedescoperi noi si noi formule de calcul numeric, chiar dac a n esent a (matematic)acestea difer a foarte put in.

Aceast a carte prezint a o viziune detaliat a asupra teoriei si practicii metodelor numerice, rod al unei activit at i de aproape 20 de ani n acest domeniu. Algorit-mii cei mai utilizat i sunt prezentat i integral. O serie de algoritmi avansat i, delarg a aplicabilitate sunt de asemenea inclusi. Autorii au ncercat o prezentare

intuitiv a a teoriei ce st a la baza metodelor numerice considerate, urm arindu-semai mult usurint a nt elegerii materialului. Locul demonstrat iilor riguroase dealtfel dicile si nu ntotdeauna eciente didactic e luat, ntr-o serie de cazuri,de observat ii crit ice si de bun simt . O observat ie de bun simt este si aceea de a face apel la mai mult a teorie atunci cand modalit at ile cunoscute au fost epuizate. Ca atare, se vor reg asi n carte si o serie de cunostint e mai avansatenecesare dezvolt arii unor metode numerice performante.

Sunt incluse capitole privind: aproximarea funct iilor, derivarea si integrarea numeric a, problemele algebrei liniare, ecuat ii si sisteme de ecuat ii neliniare, op-timizare, ecuat ii diferent iale. In prezentarea algoritmilor s-a optat pentru folo-sirea unui meta-limbaj, asem an ator celui din programul Matlab. Cititorul poatetranscrie un algoritm n limbajul de programare preferat cu usurint a. Pentru a prentimpina cererile unor utilizatori ce doresc programe surs a sau direct execu-

tabile, cartea este suplimentat a de un bogat material oferit pe Internet la adresa http://www.propulsion.pub.ro . La acest sit se pot reg asi implement ari n Pascal, FORTRAN si C++ ale celor mai utilizat i algoritmi, exemple extinse,leg aturi la alte situri de pe Internet de interes pentru analiza numeric a. Cei cu acces la Internet pot benecia de programele de instruire asistat a de calculator ce sunt disponibile la acest sit, unde este disponibil a o versiune electronic a a


4/313

acestei c art i, o serie de lucr ari de laborator si numeroase aplicat ii mici ce pot rulate direct din browser-ul preferat.

Pe tot parcursul prezent arii, elementele teoretice sunt completate cu nume-roase exemple detaliat rezolvate. Acestea provin din cele mai variate domenii:ingineria mecanic a, ingineria electric a, zic a si chimie. S-a ncercat formularea unor exemple init iale simple, ce s a se concentreze pe aspectele strict numerice,iar apoi, a unor exemple apropriate problemelor reale. Se sper a ca aceast a mo-dalitate de prezentare s a e util a at at studentului cat si practicianului metodelor numerice.

1997 Autorii


5/313

CUPRINS V

Cuprins

1 Aproximarea funct iilor de o variabil a 11.1 Aproximarea prin interpolare . . . . . . . . . . . . . . . . . . . . 2

1.1.1 Interpolarea polinomial a global a . . . . . . . . . . . . . . 31.1.2 Interpolare cu functii spline . . . . . . . . . . . . . . . . . 101.1.3 Interpolare cu functii trigonometrice . . . . . . . . . . . . 161.1.4 Interpolare n planul complex . . . . . . . . . . . . . . . . 23

1.2 Aproximarea mini-max . . . . . . . . . . . . . . . . . . . . . . . . 301.2.1 Polinoamele Cebasev . . . . . . . . . . . . . . . . . . . . . 301.2.2 Minimizarea erorii la interpolarea polinomial a . . . . . . . 321.2.3 Aproximarea aproape mini-max a unei funct ii . . . . . . . 34

1.3 Aproximarea n sensul celor mai mici p atrate . . . . . . . . . . . 361.4 Elemente de teoria aproxim arii . . . . . . . . . . . . . . . . . . . 40

1.4.1 Spatii vectoriale . . . . . . . . . . . . . . . . . . . . . . . 411.4.2 Produsul scalar si ortogonalitate . . . . . . . . . . . . . . 421.4.3 Norme, operatori si funct ionale . . . . . . . . . . . . . . . 471.4.4 Problema generala a celei mai bune aproximari . . . . . . 49

2 Derivarea si integrarea numerica 532.1 Derivarea numerica . . . . . . . . . . . . . . . . . . . . . . . . . . 53

2.1.1 Derivate folosind polinoame de interpolare . . . . . . . . . 542.1.2 Formularea operatorial a . . . . . . . . . . . . . . . . . . . 572.1.3 Polinoame de interpolare n funct ie si derivata . . . . . . 592.1.4 Derivate folosind funct ii spline . . . . . . . . . . . . . . . 612.1.5 Derivate folosind diverse aproximat ii . . . . . . . . . . . . 61

2.2 Integrarea numerica . . . . . . . . . . . . . . . . . . . . . . . . . 622.2.1 Formule Newton-Cotes nchise . . . . . . . . . . . . . . . 632.2.2 Formule de integrare deschise . . . . . . . . . . . . . . . . 682.2.3 Tehnici de atingere a unei precizii impuse . . . . . . . . . 76

3 Rezolvarea ecuat iilor neliniare 853.1 Metoda njum at atirii intervalelor . . . . . . . . . . . . . . . . . . 863.2 Procedee iterative . . . . . . . . . . . . . . . . . . . . . . . . . . 87

3.2.1 Iterat ia simpl a . . . . . . . . . . . . . . . . . . . . . . . . 873.2.2 Metoda Newton-Raphson . . . . . . . . . . . . . . . . . . 89


6/313

VI CUPRINS

3.2.3 Metoda secantei . . . . . . . . . . . . . . . . . . . . . . . 913.2.4 Metoda parabolelor tangente . . . . . . . . . . . . . . . . 93

3.3 Determinarea rad acinilor polinoamelor . . . . . . . . . . . . . . . 963.3.1 Metoda Lobacevschi-Graeffe . . . . . . . . . . . . . . . . . 963.3.2 Metode de factorizare a polinoamelor . . . . . . . . . . . 101

4 Erorile de calcul numeric 1094.1 Surse de erori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1094.2 Propagarea erorilor n calcule . . . . . . . . . . . . . . . . . . . . 111

5 Rezolvarea sistemelor liniare 1155.1 Metode directe . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

5.1.1 Metoda elimin arii a lui Gauss . . . . . . . . . . . . . . . 1165.1.2 Metoda Gauss-Jordan . . . . . . . . . . . . . . . . . . . . 1215.1.3 Propagarea erorilor la metodele de eliminare. Ranarea

solutiei . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1215.1.4 Interpretarea matriceal a a metodelor de eliminare . . . . 1245.1.5 Calculul matricei inverse . . . . . . . . . . . . . . . . . . 1255.1.6 Relat ia Sherman-Morisson . . . . . . . . . . . . . . . . . . 1285.1.7 Ranarea matricei inverse . . . . . . . . . . . . . . . . . 1295.1.8 Efectele erorilor din datele initiale . . . . . . . . . . . . . 1315.1.9 Factorizarea L U . . . . . . . . . . . . . . . . . . . . . . 1325.1.10 Descompunerea SV D . . . . . . . . . . . . . . . . . . . . 1345.1.11 Sisteme cu matrice rare . . . . . . . . . . . . . . . . . . . 135

5.2 Metode iterative . . . . . . . . . . . . . . . . . . . . . . . . . . . 1385.2.1 Metoda iterativa Jacobi . . . . . . . . . . . . . . . . . . . 1385.2.2 Metoda iterativa Gauss-Seidel . . . . . . . . . . . . . . . . 1405.2.3 Accelerarea convergentei metodelor iterative . . . . . . . . 143

5.3 Comparat ii ntre metode . . . . . . . . . . . . . . . . . . . . . . . 1455.4 Elemente de calcul matriceal . . . . . . . . . . . . . . . . . . . . 146

6 Vectori si valori proprii 1516.1 Elemente introductive . . . . . . . . . . . . . . . . . . . . . . . . 1516.2 Metode pentru cateva valori proprii . . . . . . . . . . . . . . . . . 152

6.2.1 Metoda puterii directe . . . . . . . . . . . . . . . . . . . . 1526.2.2 Metoda puterii inverse . . . . . . . . . . . . . . . . . . . . 1556.2.3 Metoda deplasarii . . . . . . . . . . . . . . . . . . . . . . 156

6.3 Determinarea tuturor valorilor si vectorilor proprii . . . . . . . . 1566.4 Metoda Danilevschi . . . . . . . . . . . . . . . . . . . . . . . . . 1576.5 Metodele QR si LR . . . . . . . . . . . . . . . . . . . . . . . . . 162

6.5.1 Rezultate teoretice preliminarii . . . . . . . . . . . . . . . 1636.5.2 Algoritmi auxiliari . . . . . . . . . . . . . . . . . . . . . . 1696.5.3 Formularea metodelor QR si LR . . . . . . . . . . . . . . 1716.5.4 Reducerea numarului de operatii la factorizare . . . . . . 1726.5.5 Accelerarea metodelor QR si LR . . . . . . . . . . . . . . 1756.5.6 Calculul vectorilor proprii . . . . . . . . . . . . . . . . . . 176


7/313

CUPRINS VII

7 Metode de optimizare 1817.1 Minimizarea n lungul unei direct ii . . . . . . . . . . . . . . . . . 1837.2 Metode de minimizare f ar a calculul derivatelor . . . . . . . . . . 1877.3 Metoda gradientului . . . . . . . . . . . . . . . . . . . . . . . . . 1907.4 Metoda Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . 1937.5 Metode cvasi-Newton . . . . . . . . . . . . . . . . . . . . . . . . . 1967.6 Metode de gradient conjugat . . . . . . . . . . . . . . . . . . . . 198

7.6.1 Rezolvarea sistemelor de ecuatii liniare folosind metodede optimizare . . . . . . . . . . . . . . . . . . . . . . . . 200

7.7 Metode specice de optimizare . . . . . . . . . . . . . . . . . . . 2047.8 Probleme de optimizare cu restrict ii . . . . . . . . . . . . . . . . 205

8 Rezolvarea sistemelor neliniare 2138.1 Iterat ia simpl a . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2148.2 Metoda iterativa Newton . . . . . . . . . . . . . . . . . . . . . . 2168.3 Metode cvasi-Newton . . . . . . . . . . . . . . . . . . . . . . . . . 2198.4 Metoda gradientului . . . . . . . . . . . . . . . . . . . . . . . . . 2218.5 Metoda hibrida . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223

9 Rezolvarea ecuat iilor diferent iale 2299.1 Considerat ii generale . . . . . . . . . . . . . . . . . . . . . . . . 2299.2 Metode cu pasi separati . . . . . . . . . . . . . . . . . . . . . . . 230

9.2.1 Formule Euler . . . . . . . . . . . . . . . . . . . . . . . . 2309.2.2 Formule Runge-Kutta . . . . . . . . . . . . . . . . . . . . 2339.2.3 Formule Runge-Kutta-Gill . . . . . . . . . . . . . . . . . . 2349.2.4 Alegerea pasului la rezolvarea ecuatiei diferent iale . . . . 235

9.3 Extrapolare Richardson. Metoda Bulirsch-Stoer . . . . . . . . . . 2389.4 Metode cu pasi legat i . . . . . . . . . . . . . . . . . . . . . . . . . 239

9.4.1 Formule explicite . . . . . . . . . . . . . . . . . . . . . . . 2399.4.2 Formule implicite . . . . . . . . . . . . . . . . . . . . . . . 241

9.5 Propagarea erorilor. Stabilitate. . . . . . . . . . . . . . . . . . . . 2439.6 Sisteme de ecuat ii diferent iale. Ecuatii de ordin superior . . . . . 247

9.6.1 Probleme cu valori init iale . . . . . . . . . . . . . . . . . . 2489.6.2 Probleme cu valori la limite . . . . . . . . . . . . . . . . . 2499.6.3 Ecuat ii diferent iale de ordin superior . . . . . . . . . . . 254

9.7 Sisteme cu scari disparate . . . . . . . . . . . . . . . . . . . . . . 255

10 Ecuat ii diferent iale cu derivate part iale 26310.1 Ecuat ii cu derivate part iale de ordinul I . . . . . . . . . . . . . . 26410.2 Ecuat ii cu derivate part iale de ordinul II . . . . . . . . . . . . . . 270

10.2.1 Ecuat ii cu derivate part iale de tip parabolic . . . . . . . . 27310.2.2 Ecuat ii cu derivate part iale de tip eliptic . . . . . . . . . . 28710.2.3 Ecuat ii cu derivate part iale de tip hiperbolic . . . . . . . 29510.2.4 Metoda caracteristicilor . . . . . . . . . . . . . . . . . . . 29710.2.5 Scheme cu diferent e nite . . . . . . . . . . . . . . . . . . 300


8/313

1

Capitolul 1

Aproximarea funct iilor de ovariabila

Problema aproxim arii unei funct ii de o variabil a se poate pune n situat iidiverse, urm atoarele doua ind mai frecvente:

1. functia este cunoscuta, dar are o form a complicat a, dicil de manipulat ncalcule (spre exemplu pentru operat ii de derivare, integrare, etc.);

2. functia nu este complet cunoscut a, ind date numai valorile ei pe o multimediscret a si nit a de puncte.

In primul caz, aproximarea se poate face, n principiu, oric at de exact, res-trict iile ind legate de conditia ca funct ia care aproximeaza sa e cat mai simpl a.In al doilea caz informat iile sunt reduse si se completeaz a cu presupuneri supli-mentare, privind gradul de regularitate al funct iei (continuitatea funct iei si aderivatelor sale, etc.). In ambele cazuri, este important a alegerea unui criteriu de aproximare .

Fie [a, b] R un interval pe dreapta real a si x i , i {1, 2, . . . , N }, N N ,o retea de puncte de diviziune ale acestui interval, xi [a, b], x1 = a, x N = b.Punctele de diviziune se numesc noduri . Presupunem date valorile n noduri ale


9/313

2 1. Aproximarea functiilor de o variabil a

unei funct ii reale f

yi = f (x i ), i 1, N . (1.1)

Not am cu g(x) funct ia cu care vrem s a aproxim am pe f (x) pe intervalul dat.Iat a cateva criterii de aproximare:

a) Interpolare . In acest caz, funct ia mai simpla g(x) este determinat a dincondit ia ca s a ia aceleasi valori n noduri

g(xi ) = yi , i 1, N . (1.2)

Criteriul de aproximare prin interpolare presupune tacit c a nodurile ( x i , yi )sunt cunoscute exact. Dac a din diverse motive cel mai adesea datorit aunui procedeu de masurare nodurile sunt afectate de erori atunci criteriulde interpolare este inadecvat.

b) Minimizarea abaterii maxime. Se impune condit ia ca abaterea maxim a sae minima pe intervalul ales, adica

maxx [a,b ] |f (x) g(x)| = minim. (1.3)

Relat ia ( 1.3) are analogul discret

maxi 1,N |yi g(x i )| = minim. (1.4)

Aproximarea facut a pe baza criteriului de mai sus se numeste aproximaremini-max .

c) Minimizarea sumei p atratelor abaterilor n noduri . In acest caz se impuneca

S =n

i=1

(yi g(x i ))2 = minim . (1.5)

Se observa ca, n cazul interpol arii, aceast a suma este chiar nula, adic a arecea mai mica valoare posibil a. Totusi, aceast a observat ie nu face superuucriteriul ( 1.5) care este mai general si permite tratarea datelor cunoscuteincert, asa cum se va vedea mai departe. Metoda care foloseste acestcriteriu este nt alnit a sub numele de metoda celor mai mici p atrate .

1.1 Aproximarea prin interpolare

Presupunand c a nodurile x i sunt distincte , condit ia de interpolare ( 1.1) repre-zint a un sistem de N condit ii si va duce n general la un sistem de N ecuat iicu N necunoscute. Considerentele de simplitate amintite mai sus ne sugereaz a


10/313

1.1. Aproximarea prin interpolare 3

c-alegerea formei funct iei de aproximare sa e facut a astfel nc at sistemul decondit ii sa conduc a la ecuat ii liniare. O posibil a alegere este urmatoarea: seia un set de N funct ii simple, cunoscute, gk (x), k {1, 2,...,N }si un set deN parametrii nedeterminat i (scalari) ak , k {1, 2,...,N }, n funct ie de care sescrie aproximanta g(x)

g(x) =N

k=1

ak gk (x) . (1.6)

Deoarece nu ne-am nscris ntr-un formalism riguros, vom face unele observat iide bun simt . Astfel, am ales N parametri nedeterminat i ak , deoarece avemN condit ii. Pe de alta parte, setul de funct ii gk (x) trebuie s a contin a elementedistincte, astfel nc at introduse n forma ( 1.6), num arul de parametri s a nu sereduc a1 . Intr-un limbaj mai riguros, se spune c a cele N funct ii cunoscute gktrebuie s a e liniar independente . In lipsa altei analize, ne putem limita lafunct ii despre care stim ca au aceasta proprietate. Un astfel de set l reprezint amonoamele xk 1 , k {1, 2,...,N }, n care caz funct ia de interpolare este unpolinom de gradul N 1

g(x) =N

k=1

ak xk 1 . (1.7)

Alte seturi de functii pot funct iile trigonometrice, exponent iale, etc., pe carele vom trata ulterior.

1.1.1 Interpolarea polinomiala global a

Revenind la forma ( 1.7) se pune problema gasirii coecient ilor ak din condit ia

de interpolare, adica a rezolvarii sistemului de ecuatii liniareN

k =1

ak xk 1 = yi , i 1, N . (1.8)

Daca N este mare, rezolvarea sistemului ( 1.8) este dicil a sau cel put in necon-venabil a. In orice caz, nodurile x i ind distincte, sistemul de ecuat ii (1.8) esteun sistem cu determinant Vandermonde diferit de zero si are o solut ie unic a,bine determinat a pentru coecientii ak . In consecint a, oricare ar calea pe carese construieste efectiv polinomul de interpolare ( 1.7), acesta este unic pentru ofunct ie si o diviziune data. Aproximarea efectuat a este o interpolare global a nsensul ca se foloseste un singur polinom pe tot intervalul [ a, b].

Forma Newton a polinomului de interpolare.O modalitate convenabil a de constructie a polinomului de interpolare l consti-tuie polinomul lui Newton cu diferent e divizate . Fiind date perechile de puncte

1 Spre exemplu, dac a g2 = g 1 , = 0 atunci a 1 g1 + a 2 g2 = ( a 1 + a 2 )g1 , deci n loc de doiparametri independent i a 1 si a 2 ar apare doar o combinat ie a 1 = a 1 + a 2 .


11/313


(x i , yi ), se introduc urmatoarele rapoarte denumite diferent e divizate (DD)

DD( x2 , x1) = y2 y1x2 x1,

DD( x3 , x2 , x1) =DD( x3 , x2) DD( x2 , x1)

x3 x1,

DD( xN , xN 1 , . . . , x 1) =DD( xN , . . . , x 2) DD( xN 1 , . . . , x 1)

xN x1. (1.9)

Diferent ele divizate care se construiesc folosind k +1 puncte se numesc diferent edivizate de ordinul k. Se poate demonstra prin induct ie matematic a, urm atoareaexpresie a diferent ei divizate de ordinul N 1:

DD( xN , xN 1 , . . . , x 1) =N

i=1

yi N

j =1

(x i x j ) . (1.10)

Semnul denot a omiterea factorului j = i din produs. Relatia ( 1.10) indsimetric a, rezult a ca valoarea diferentei divizate nu depinde de ordinea n careluam punctele x i . Din punct de vedere practic, este mai comod ca diferent eledivizate s a nu se calculeze cu formula ( 1.10) ci recursiv, alc atuindu-se un tabel(Tab. 1.1).

Funct ia de aproximare g(x) este un polinom de gradul N 1, pe care l vomnota cu pN 1(x). Vom scrie

f (x) = pN 1(x) + RN 1(x) , (1.11)

unde RN 1(x) este restul sau eroarea de aproximare la interpolarea polinomial a.Pe de alt a parte, t in and seama de denitia diferent elor divizate, se poate scrie

f (x) = y1 + ( x x1)DD( x, x 1) =y1 + ( x x1)DD( x2 , x1) + ( x x1)(x x2)DD( x, x 2 , x1) .(1.12)

In relat ia ( 1.12) ultima diferenta divizat a este formata cu punctul curent x.Continuand procedeul ( 1.12) pan a se iau n consideratie toate nodurile, rezult a

pN 1(x) = y1 + ( x

x1)DD( x2 , x1) + ( x

x1)(x

x2)DD( x, x 2 , x1) +

(1.13). . . + ( x x1)(x x2) . . . (x xN 1)DD( xN , xN 1 , . . . , x 1) ,

RN 1 =N

i=1

(x x i ) DD( x, x N , xN 1 , . . . , x 1) . (1.14)


12/313


Se verica direct din ( 1.14) ca restul se anuleaza n noduri ( RN 1(x i ) =0, i = 1 , 2,...,N ) si deci pN 1(x) dat de ( 1.13) este polinomul de interpolarecorespunz ator celor N puncte date. Forma ( 1.13) se numeste polinomul lui Newton cu diferent e divizate.

Pentru scrierea polinomului ( 1.13) se alcatuieste mai int ai tabloul diferentelordivizate de diverse ordine. Nu este necesar ca tabloul s a cuprind a toate diferenteledivizate posibile. Este sucient s a proced am n ordinea nodurilor, din aproapen aproape, conform denit iilor (1.9). Un exemplu de calcul este prezentat ntabelul 1.1 unde s-a notat cu DD i diferent a divizata de ordinul i (i = 1 , 2, 3).Polinomul obt inut este

p3(x) = 2 + ( x 1) (1) + ( x 1)(x 2) (2) + ( x 1)(x 2)(x 3) (1) .In ceea ce priveste restul RN 1(x) (eroarea la interpolare), se poate face o

evaluare daca avem informat ii suplimentare referitoare la funct ia aproximat af (x) si la derivatele sale. In acest scop, consideram funct ia ajutatoare Q(t)denit a prin relat ia

Q(t) = f (t) pN 1(t) N

i =1(t x i ) DD( x, x N , xN 1 , . . . , x 1) . (1.15)

Se observa ca funct ia Q(t) se anuleaz a pentru t = x si t = x i , i = 1 , 2,...N , adicaare N + 1 zerouri. Presupun and c a f (t) este derivabila de un num ar convenabilde ori, putem aplica functiei Q(t) si derivatelor sale teorema lui Rolle. Rezult a,succesiv, c a derivata Q (t) are cel put in N zerouri, derivata Q (t) are cel put inN 1 zerouri s.a.m.d., astfel ca derivata de ordinul N are cel put in un zero peintervalul ( a, b). Fie t = acest zero. Deriv and relat ia ( 1.15) de N ori, si facandt = , se obtine

f (N ) () = N ! DD( x, x N , xN 1 , . . . , x 1) ,relat ie din care putem deduce expresia diferent ei divizate de ordinul N n funct iede derivata de ordinul N . In acest fel, restul la interpolare ( 1.14) cap at a forma

RN 1(x) =N

i =1

(x x i ) f (N ) ()/N ! . (1.16)

Prezent a produselor (x x i ), sugereaz a ca restul este mai mic (n modul)cand punctul curent x este centrat pe intervalul care cont ine diviziunea si maimare c and x este luat spre marginile intervalului sau n afara lui acest ultim

caz este denumit extrapolare . Deoarece derivata de ordinul N a functiei npunctul nu este accesibil a (din diverse motive), evalu ari ale restului se potface presupunand c a, la schimbarea diviziunii, punctul (necunoscut) nu sedeplaseaz a mult, astfel ncat derivata respectiv a sa e aproximativ constant a,si refacand calculul pentru o nou a diviziune a intervalului se poate testa astfelsi sensibilitatea erorii la includerea de noi puncte de interpolare.


13/313


Tabelul 1.1: Un tabel de diferent e divizate

x i yi DD1 DD2 DD3 DD41 2 2 3 1 3 0 -3 -2 5 6 3 2 1 4 4 2 1 -3/2 -5/6

De regul a, este nerecomandabil a utilizarea unui num ar mare de noduri lainterpolare, mai ales daca se intent ioneaz a calcularea unor derivate cu ajutorulacestui polinom. Uzual, nu se folosesc toate cele N noduri, ci doar 3-5 noduricele mai apropriate de punctul n care se cere valoarea funct iei. In consecint a,exist a posibilitatea unor variante chiar c and nu putem ndesi ret eaua prin ale-

gerea altor noduri.Interpolarea polinomial a apare deseori ca o component a a altor algoritminumerici cum ar integrarea sau derivarea numeric a a unor funct ii. In acesteaplicat ii se consider a deseori cazul diviziunilor egale

x i +1 x i = h, i = 1 , 2, . . . , N 1 ,h ind pasul ret elei . Se introduc operatorii si denumit i diferent a la dreapta

respectiv, diferent a la st anga , prin relat iile

f (x) = f (x + h) f (x) , (1.17)f (x) = f (x) f (x h) . (1.18)

Rezultatul aplicarii operatorului sau asupra lui f (x) se numeste diferent a nit a (de ordinul I). Pentru n ntreg, se deneste un operator de translat ie, E ,prin relat ia

E n f (x) = f (x + nh ), n Z . (1.19)

Avem E 1f (x) = f (x + h), E 0f (x) = f (x), E 1f (x) = f (x h). Se observ a cantre operatorii , si E exist a relat iile = E E 0 , = E 0 E 1 . (1.20)

Diferent ele divizate se pot exprima n funct ie de diferent ele nite si de pasul h

DD( x2 , x1) = [f (x1 + h) f (x1)] /h = [ f (x1)] /h , (1.21)

DD( xN , xN 1) = [f (xN ) f (xN h)] /h = [ f (xN )] /h . (1.22)Prin inductie, se demonstreaza usor c a

DD( xN , xN 1 , . . . , x 1) = N 1f (x1)

(N 1)!hN 1=

N 1f (xN )(N 1)!hN 1

, (1.23)


14/313


unde exponentul indic a aplicarea repetat a a operatorului.Pun and variabila curent a x sub forma

x = x i + h, [0, N 1] , (1.24)se poate obt ine expresia polinomul de interpolare Newton cu diferent e nite la dreapta

pN 1(x) = y1 + y1 + 12 ( 1) 2y1 + . . . + C N 1 N 1y1 , (1.25)unde C k , k = 0 , 1, . . . , N 1 sunt coecient ii binomiali. Restul RN 1(x) cap at aforma

RN 1(x1 + h ) = hN C N f (N ) () . (1.26)

Calculul se face alc atuind un tablou al diferent elor nite, similar cu tabloul

diferent elor divizate.In mod asem an ator, notand

x = xN + h, [N + 1 , 0] , (1.27)se obtin expresii cu diferent e la st anga

pN 1(x) = yN + yN + 12 ( + 1)2yN + . . . + ( 1)N 1C N 1 N 1yN

RN 1(xN + h ) = ( 1)N hN C N f (N ) () . (1.28)Forma Lagrange a polinomului de interpolare.

Polinomul de interpolare Lagrange se scrie alegand funct iile gk (x) din relat ia(1.6) sub forma unor polinoame denumite polinoame Lagrange si notate cu

Lk (x), k = 1 , 2,...,N . Aceste polinoame au expresiile

Lk (x) =N

j =1

x xjxk x j

, k 1, N , (1.29)

unde produsul se ia pentru j = k. Se observa direct din ( 1.29) ca

Lk (x j ) = 0 dac a x j = xk ; Lk (x j ) = 1 dac a xj = xk . (1.30)Polinomul de interpolare Lagrange se scrie

pN 1(x) =N

k=1ykLk (x) , (1.31)

deci coecientii ak din expresia ( 1.6) sunt chiar valorile functiei f (x) n noduri,ak = yk . Se verica direct, t in and cont de propriet at ile (1.30) ca pN 1(x i ) =yi , i = 1 , 2,...,N .


15/313


16/313


Figura 1.1: Polinoamele de interpolare pN 1 (x) cu N = 6 , 11 suprapuse peste f (x).

Aspecte practice ale interpol arii polinomiale.

Evident, polinoamele Newton si Lagrange difer a numai prin forma, pentruaceeasi retea restul ind acelasi n ambele cazuri. Din punct de vedere al calcu-lului numeric, este preferat a folosirea polinomului Newton ce necesita un num arde operat ii aritmetice mai mic, de O(3N 2 / 2) fat a de O(4N 2) pentru polinomulLagrange. Necesarul de memorie este acelasi pentru ambii algoritmi. Pentrupolinomul Newton ar parea c a este nevoie de o matrice suplimentar a pentrutabelul de diferente divizate. Insa din tabelul de diferente divizate se folosesc

efectiv doar N coecienti exist and posibilitatea refolosirii celulelor de memorieprecum n algoritmul 2

d y j = 2 : N k = N : 1 : j[dk (dk dk 1)/ (xk xk j 1)

(1.32)

n urma caruia diferent ele divizate de pe diagonala tabelului se obt in n vectoruld ce a fost initializat cu ordonatele y. Aceasta este partea cea mai costisitoarea interpol arii Newton necesitand O(3N 2 / 2) operat ii aritmetice. Evaluarea po-linomului ntr-un punct u se face ecient prin schema lui Horner

S dN j = ( N 1) : 1 : 1[S dj + ( u xj ) S

(1.33)

2 Am folosit notat ia Matlab j = j ini : pas : j fin pentru bucle: variabila j este initializat a lavaloarea j ini si apoi incrementat a cu pas . Instruct iunile din bucl a, delimitate de [ se execut arepetat p an a c and j > j fin Dac a pas nu este precizat, precumn j = j ini : j fin , se subnt elege pas = 1.


17/313


Figura 1.2: Aproximarea funct iei lui Runge f (x) printr-o linie poligonal a s(x).

necesit and numai O(3N ) operat ii. Evaluarea polinomului Lagrange ntr-unpunct u necesit a O(4N 2) prin

S 0k = 1 : N P 1 j = 1 : k 1[P P (u x j )/ (xk x j )

j = k + 1 : N [P P (u x j )/ (xk x j )S S + yk P

1.1.2 Interpolare cu funct ii spline

Am vazut c a interpolarea polinomial a global a , pe tot intervalul [ a, b], nu convergentotdeauna. Desigur, dac a am reusi s a micsor am diviziunea f ar a a modica gra-dul polinomului de interpolare, rezultatul ar putea modicat. Spre exemplu,aproximarea unei funct ii derivabile cu o linie poligonal a se poate face oric at debine cand numarul de laturi ale poligonului creste innit (g. 1.2). Evident,funct ia poligonal a nu se identic a cu o funct ie de gradul 1 deoarece depindesi de diviziunea aleasa. Acest exemplu conduce la ideea de interpolare polino-miala pe port iuni, la care pe ecare subdiviziune a intervalului [ a, b] denim unalt polinom de interpolare. Funct ia poligonal a este unul dintre exemplele celemai simple ale acestui tip de interpolare prin funct ii spline 3 . Aceste funct ii

sunt caracterizate prin formele lor pe subintervalele dintre dou a noduri (carepot diverse funct ii cunoscute) si prin anumite condit ii de racordare n noduri.In cele ce urmeaz a, vom considera doar cazul funct iilor spline polinomiale f ar a decient a .

3 Se citeste splain.


18/313


Figura 1.3: Subintervalele de denire a unei funct ii spline.

Denit ie. Fie [a, b] R un interval pe dreapta real a si x i , i = 1 , 2,...,N oretea de puncte de diviziune ( x1 = a, x N = b). Not am cu I i subintervalele[x i , x i +1 ). Funct ia s : [a, b] R se numeste funct ie spline polinomial a deordinul m daca

1. restrict iile ei pe subintervalele I i sunt polinoame de gradul m, s|I i = pm,i ;2. s este derivabil a de m 1 ori pe intervalul [ a, b], s C (m 1) [a, b].

A doua condit ie cont ine n sine conditia de racordare n noduri

p(k )m,i (x i+1 ) = p(k )m,i +1 (xi +1 ), k = 0 , 1, . . . , m 1 , (1.34)adica la frontiera x i +1 dintre doua subintervale, polinomul din st anga pm,i siprimele sale m 1 derivate trebuie sa aibe aceleasi valori cu ale polinomuluidin dreapta, pm,i +1 . In afara intervalului [ a, b] functia s se poate prelungi prinpolinoame de grad m. Condit iile de racordare n noduri pot sl abite, astfelnc at funct ia s sa nu mai e de clas a C (m 1) pe tot intervalul [ a, b], ci sa ederivabil a de mai put ine ori pe diverse subintervale. In acest caz, obt inem funct ii spline cu decient a .

Funct ia spline de ordinul nt ai (linia poligonala).

Funct ia spline este formata din segmente de ecuatie

p1,i (x) = yi + m i (x x i ), x [x i , x i +1 ), (1.35)m i = ( yi+1 yi )/h i , hi x i +1 xi , (1.36)

m i reprezentand panta pe intervalul I i (vezi g. 1.3). Funct ia spline de ordinulntai este simpl a, dar nu furnizeaza derivata functiei interpolate.


19/313


Funct ia spline de ordinul doi.

Funct ia este formata din segmente de parabol a, racordate n noduri p ana laderivata de ordinul 1, inclusiv

p2,i (x) = yi + m i (x x i ) + a i (x x i )2 , x [x i , x i +1 ), i 1, N 1 . (1.37)Forma polinomiala (1.37) satisface automat condit ia p2,i (x i ) = yi prin modulde scriere. Condit iile de racordare

p2,i (x i +1 ) = yi +1 (1.38)

conduc la urmatoarele ecuatii pentru coecientii a i

a i = ( yi +1 yi )/h 2i m i /h i , i 1, N 1 . (1.39)Eliminarea lui a i din condit iile de racordare

p2,i (x i +1 ) = yi +1 , p2,i (x i +1 ) = p2,i +1 (x i +1 ) ,

(1.40)

care se pot scrie n nodurile x i , i = 1 , 2, 3, . . . , N 2 conduce la sistemulm i + m i +1 = 2( yi +1 yi )/h i , i 2, N 1 , (1.41)

ce trebuie completat cu o singur a condit ie. Spre exemplu, se poate da panta launul din capetele intervalului ( m1 sau mN ). Necesitatea conditiei suplimentareprovine din faptul ca nu mai putem impune condit ia de racordare n derivat a nnodul xN . In ambele cazuri, sistemul devine determinat:

1. m1 dat duce la substituirea

m i +1 = 2( yi +1 yi )/h i m i , i 1, N 1 ; (1.42)2. mN dat permite retrosubstituirea

m i = 2( yi +1 yi )/h i m i+1 , i N 1, 1 . (1.43)Funct ia spline de ordinul trei sau cubic a.

Este una din cele mai utilizate funct ii spline, av and derivate continue p ana laordinul doi inclusiv, ceea ce permite calculul razei de curbur a. Din condit iile decontinuitate pe noduri p ana la derivata de ordinul 2 inclusiv

p3,i (x i +1 ) = yi +1 , p3,i (x i +1 ) = p3,i +1 (x i +1 ) , p3,i (x i +1 ) = p3,i +1 (x i +1 ) ,

(1.44)


20/313


pentru i = 1 , 2, 3, . . . , N 2, se deduc coecient ii polinomului de gradul 3 p3,i (x) = yi + m i (x x i ) + bi (x x i )

2+ a i (x xi )

3(1.45)

care reprezinta comportarea funct iei spline pe ecare subinterval ( x i , x i+1 ), i1, N 1

a i = ( m i +1 + m i )/h 2i 2(yi +1 yi )/h 3i , (1.46)

bi = 3( yi +1 yi )/h 2i (m i +1 + 2 m i )/h i . (1.47)Pantele pe noduri, m i , sunt date n acest caz de sistemul

i m i 1 + 2 m i + i m i = di , i 2, N 2 (1.48)ce trebuie completat cu dou a condit ii, pentru a suplini condit iile de racordaren prima si a doua derivat a ce nu mai pot scrise n xN . S-au facut notatiile

i h i / (h i 1 + h i ), i 1 i , hi x i +1 x i , (1.49)

di 3 [ i (yi +1 yi )/h i + i (yi yi 1)/h i 1] .Sistemul de condit ii de racordare impuse lasa de data aceasta dou a grade de

libertate, ce pot precizarea pantelor la capete, m1 si mN , sau, mai general,precizarea unei relatii, n general liniar a, ale acestor pante cu pantele vecine, deforma

2m1 + 1m2 = d1N mN 1 + 2 mN = dN .

. (1.50)

In relat iile (1.50), coecientii 1 , d1 , N , dN sunt dati prin natura condit iilorpuse la capetele intervalului, deci nu sunt dedusi din relat iile (1.49) care nu sunt denite pentru i = 1 si i = N . Spre exemplu, a da pantele m1 si mN revine la a impune

1 = 0 , d1 = 2 m1 , N = 0 , dN = 2 mN .

Sistemul de N ecuat ii cu N necunoscute Am = d, obtinut prin reuniunea egali-t at ilor ( 1.48) si (1.50), are matrice tridiagonal a. Ne punem mai nt ai problemaexistent ei unei solut ii. Elementele din matricea A rezultate din condit iile deracordare ( 1.48) sunt diagonal dominante pe linii, adica 2 > |i | + | i | = 1.Cum o matrice diagonal dominant a este, n general, inversabil a, este sucientca si condit iile suplimentare ( 1.50) sa pastreze aceasta proprietate. Practic,problema este rezolvabila dac a impunem conditii necontradictorii si distincte.

Forma tridiagonal a a matricei A permite o rezolvare foarte ecient a prindescompunerea matricei ntr-un produs de dou a matrice bidiagonale 4 A = L R

4 Un caz particular al factoriz arii Doolittle ce va prezentat a n capitolul 5, cunoscut caalgoritmul lui Thomas .


21/313


sau explicit

A =

1 0 0 0l2 1 0 00 l3

. . ....

......

.... . . 1 0

0 0 lN 1

r 1 1 0 00 r 2

. . . 0 0

0 0. . . N 2 0

......

... r N 1 N 10 0 0 r N

.

Coecient ii necunoscut i r i , i = 1 , 2,...,N si li , i = 2 , 3,...,N se determin a prinidenticarea elementelor din A si din matricea produs. Elementele supradiago-nalei din matricea R au fost deja identicate. Determinarea pantelor se bazeaz ape asociativitatea produsului de matrice

(L

R)

m = L

(R

m) = d .

Introducem vectorul z = R m. Etapele algoritmului sunt1. factorizarea A = L R

r 1 2;i = 2 : N [li i /r i 1 ; r i 2 li i 1

2. rezolvarea sistemului bidiagonal L z = dz1 d1i = 2 : N

[zi di li zi 13. rezolvarea sistemului bidiagonal R m = z

mN zN /r N i = ( N 1) : 1 : 1[m i (zi i m i +1 )/r iObservat ii. 1. Un caz particular important de funct ie spline de ordinul treieste cel al interpol arii spline cubice naturale denit prin conditiile n capete

s (x1) = s (xN ) = 0

ceea ce conduce la

2m1 + m2 = 3( y2 y1)/h 1 ,

mN 1 + 2 mN = 3( yN yN 1)/h N 1 , (1.51)


22/313


23/313


Figura 1.4: Logaritmul zecimal al erorii relative la interpolarea spline a f, f , f .

2. extrapolarea liniar a a pantelor adiacente m 1 = m 2 , mN 1 = m N ceea ceconduce la 1 = 2, d1 = 0 , N = 2, dN = 0 ;

3. condit ii naturale ( 1.52).Calitatea aproxim arii se apreciaza prin evaluarea erorilor relative

k =1000

i =1

f ( k ) (u i ) s ( k ) (u i ) / f ( k ) (u i ) , k = 0 , 1, 2

pentru diverse valori ale lui N unde {u i , i 1, 1000} este o diviziune echidistant ana a intervalului [ 5, 5]. Rezultatele sunt prezentate n gura ( 1.4). Se observ aca interpolarea este convergent a, eroarea relativ a scazand rapid cu cresterea lui N .Pentru un num ar mic de puncte toate condit iile de cap at dau rezultate comparabile. Pemasur a ce diviziunea devine sucient de n a pentru a descrie precis variat iile funct iei,condit iile exacte dau eroarea minim a dup a cum era de asteptat, urmate de condit iilenaturale si apoi de cele de extrapolare a pantei. Se remarc a pierderea a 1,2 ordine deprecizie a aproximarii n urma deriv arii.

1.1.3 Interpolare cu funct ii trigonometrice

Interpolarea polinomial a nu este adecvata aproxim arii tuturor tipurilor de funct ii.Vom considera acum o alt a alegere a funct iilor liniar independente din ( 1.6),anume funct iile trigonometrice (g. 1.5)

cos(2kx ), k 0, N ; sin(2mx ), m 1, N 1 . (1.54)Aceast a baz a este deosebit de ecace n aproximarea funct iilor periodice f (x) =f (x + 1). Functiile periodice cu o alt a perioad a f (z) = f (z + T ) pot adusela forma anterioara prin transformarea x = z/T . In cazul interpolarii cu funct ii


24/313


Figura 1.5: Primele 11 functii din baza trigonometric a.

trigonometrice, funct iile sinus si cosinus mpreun a formeaza o baza5 . Avem unnum ar de 2N functii n aceasta baz a. Ca atare, vom considera un num ar par2N de puncte de diviziune echidistante pe intervalul [0 , 1]

x j = j/ 2N, j 0, 2N 1 . (1.55)Se verica direct c a setul ( 1.54) prezint a urm atoarele propriet ati de ortogonali-tate pe multimea discreta de puncte

{x i

}=

{0, 1/ 2N, 2/ 2N, ..., (2N

1)/ 2N

}2N 1

j =0

cos2kx j cos2mx j =0, k = mN, k = m = 0 , N 2N, k = m = 0 , N

(1.56)

2N 1

j =0

sin2kx j sin2mx j =0, k = mN, k = m ;

2N 1

j =0

cos2kx j sin2mx j = 0 ,

k 0,N, m 1, N 1 . Demonstrat ia se construieste prin transformareaproduselor de functii trigonometrice n sume de sinus si cosinus. Acestea se potnlocui cu functii exponent iale, sin x = ( e

ix

e ix

)/ 2i, cosx = ( eix

+ e ix

)/ 2,rezult and progresii geometrice simplu de nsumat (vezi si 1.1.4).5 Strict vorbind doar mult imea innit a {1, cos x, sin x, cos2 x, sin2 x , . . . } formeaza o baz a a

spat iului de funct ii. Pastrarea unui num ar nit de funct ii conduce la aparit ia unor erori cevor considerate mai jos.


25/313


Utiliz and setul de functii de baz a (1.54), aproximanta prin interpolare sescrie sub forma polinomului Fourier

g(x) = 12 a0 +N 1

k=1

[ak cos(2kx ) + bk sin(2kx )] + 12 aN cos2N x , (1.57)

ce satisface 2N condit ii de interpolare

g(x j ) = f (x j ) yj , j 0, 2N 1 .Coecient ii ak , bk se determin a prin utilizarea propriet at ilor de ortogonalitate(1.56). Prezentam calculul doar pentru coecient ii bk . Polinomul ( 1.57) seevalueaz a n punctele x j , relat ia obt inuta se amplica cu sin 2mx j iar apoi secalculeaz a suma de la j = 0 la j = 2 N 1

2N 1

j =0yj sin2mx j = a

0

22N 1

j =0sin2mx j +

N 1

k=1ak

2N 1

j =0(cos2kx j sin2mx j ) +

bk2N 1

j =0(sin2kx j sin2mx j ) +

aN 2

2N 1

j =0(cos2N x j sin2mx j ) .

Se poate observa schimbarea ordinii de nsumare din relat ia de mai sus ce per-mite aplicarea ( 1.56) obtin andu-se

2N 1

j =0g(x j ) sin(2mx j ) =

2N 1

j =0yj sin(2mx j ) = bm N .

Un calcul analog pentru ceilalt i coecient i conduce la relat iile

ak =1N

2N 1

j =0

yj cos2kx j , bm =1N

2N 1

j =0

yj sin2mx j (1.58)

cu k 0,N, m 1, N 1.In aplicat ii, coecient ii ak , bk se evalueza mult mai economic dec at princalculul direct al sumelor de mai sus prin folosirea transform arii Fourier rapideprezentate n 1.1.4. Se poate lesne observa din ( 1.58) ca vom avea tot i ak = 0pentru functii impare f (x) = f (x) si tot i bm = 0 pentru funct ii pare f (x) =f (x).

Aparit ia unei oarecare asimetrii termenii n cos sunt mai numerosi dec at cein sin este legata de alegerea unui numar par de 2 N intervale n care se divide

perioada funct iei. Daca se aleg 2N + 1 intervale, forma funct iei de interpolareeste

g(x) = 12 a0 +N

k=1

[ak cos(2kx ) + bk sin(2kx )] ,


26/313


coecientii ak si bk ind dat i de

ak = 22N + 1

2N

j =0yj cos2kx j , bm = 22N + 1

2N

j =0yj sin2mx j ,

cu k 0,N, m 1,N, x j = j/ (2N + 1).

Convergent a interpol arii trigonometrice.

Se pot determina exprim ari ale erorii de interpolare n genul restului RN 1(x)de la interpolarea polinomial a si pentru interpolarea trigonometric a. Expresiileobtinute sunt ns a sub o forma integral a greu utilizabil a practic. Vom preferao discut ie mai put in formal a a erorii si convergentei procedeului. Daca funct iaf admite o dezvoltare n serie Fourier m arirea numarului de noduri conducela aproximat ii din ce n ce mai bune. Intr-adevar, ridic and expresia ( 1.57) lapatrat membru cu membru, nsum and valorile pe noduri si tin and seama derelat iile de ortogonalitate ( 1.56), se obtine

14

a20 +12

N 1

k=1

(a2k + b2k ) +

14

a2N =1

2N

2N 1

j =0

y2j (1.59)

relat ie denumita egalitatea lui Parseval discret a . Cand N creste, suma dinmembrul drept se aproprie de integrala

10 y

2dx. Daca integrala este marginit a6suma este de asemenea marginit a. Ca urmare, seria patratelor coecientilor esteconvergenta, ceea ce arat a ca a2k , b

2k devin din ce n ce mai mici c and N creste.

Interpolarea trigonometric a este asadar convergent a pentru functii f continuesau cu un numar nit de discontinuit at i acestea ind conditii suciente pentrua asigura existenta integralei anterioare. Stabilirea convergent ei este important apentru validarea interpol arii trigonometrice. In aplicat ii nsa mai apar si alteaspecte ale comportarii erorii la interpolare. Vom considera doar dou a maiimportante: rapiditatea convergent ei si efectul considerarii doar a unui numarnit de funct ii trigonometrice n dezvoltarea ( 1.57).

Far a a ncerca o denitie formal a, vom spune c a o interpolare trigonome-tric a este rapid convergent a dac a num arul de termeni N necesar realiz arii uneiprecizii impuse a aproximarii este mic. Int elesul cuv antului mic depindede aplicat ie, dar un domeniu orientativ ar 2 N 128. Urmarirea g. 1.5sugereaz a ca includerea mai multor termeni n polinomul Fourier ( 1.57) permitedescrierea unei funct ii cu variat ii mai rapide pe intervalul [0 , 1]. Fie x cea maimica distant a de pe abscis a pe care funct ia f are o variat ie semnicativa. Deo-arece f are perioada 1, spunem ca = 1 / x este frecvent a variatiilor celor mai

rapide ale funct iei. Pentru a descrie variat iile date de frecvent a cea mai rapidaa funct iei f polinomul ( 1.57) trebuie s a contin a un num ar de termeni N .Acest rezultat este cunoscut sub numele de criteriul Nyquist ce rezult a dintr-un rezultat mai general denumit teorema de esantionare Shannon . Observat i ca

6 Intr-o formulare riguroas a dac a y = f (x ) este patrat integrabil a pe [0 , 1], ceea ce se scrief L 2 [0, 1].


27/313


Figura 1.6: Comportarea polinoamelor trigonometrice cu N = 8 , 16, 32, 64 la inter-polarea functiei treapt a. Pe m asur a ce N creste, interpolarea se aproprie mai mult defunct ia exacta far a n s a a elimina oscilat iile n zona punctelor de discontinuitate.

pentru a avea N frecvent e n polinomul ( 1.57) este nevoie de 2N puncte ( x j , yj ).De aici o formulare echivalenta a criteriului Nyquist este ca num arul de noduri s a e minim de dou a ori frecvent a cea mai rapid a .

Daca f are variat ii lente atunci este mic si num arul de termeni din poli-nomul ( 1.57) este de asemenea mic. Funct iile cu variat ii rapide au ns a marenecesit and un numar mare de termeni n polinomul ( 1.57). Cea mai rapid afrecvent a de variat ie posibil a a unei funct ii ar

ceea ce corespunde

la x = 0, adic a f sa prezinte discontinuit at i. Cum nu vom putea ndeplininiciodat a criteriul Nyquist N pentru functii discontinue, interpolarea tri-gonometric a va avea erori mai mari n asemenea cazuri. Erorile ce apar suntconcentrate n jurul discontinuit atilor, comportare cunoscut a sub denumirea de fenomenul Gibbs. Un exemplu faimos al fenomenul Gibbs este aproximarea unuisemnal dreptunghiular y(x) = 1 pentru n < x < n + 1 / 2, y(x) = 1 pentrun + 1 / 2 < x < n + 1 si y(n + 1 / 2) = 0 cu n N , exemplu prezentat n g.1.6. Ne asteptam asadar ca interpolarea trigonometric a sa e lent convergentapentru functii discontinue si rapid convergent a pentru functii netede, cu variat iilente.

Sa presupunem acum ca nu este ndeplinit criteriul Nyquist 7 si am luat unnum ar prea mic N < de termeni n dezvoltarea ( 1.57). Ne punem problemadaca coecientii ak , bk k N determinati prin relat iile (1.58) sunt corect i, adic aau aceleasi valori ca n cazul n care criteriul Nyquist ar satisf acut. R aspunsuleste negativ . Sa refacem calculul anterior ce a furnizat valorile coecient ilorbk , de data aceasta pentru dezvoltarea Fourier complet a a funct iei f ce are

7 Deoarece funct ia f este n general necunoscut a si ca atare nu cunoastem frecvent a .


28/313


coecientii exact i k , k

f (x) = 12 0 +

k =1

[ k cos(2kx ) + k sin(2kx )] . (1.60)

Ca mai nainte, evalu am (1.60) n x j , nmul t im cu sin2 mx j si nsumam de la j = 0 la j = 2 N 1. Urmarim doar termenii cu produse de funct ii sinus ceilalti dau contribut ii nule conform ( 1.56)

2N 1j =0 f (x j )sin2mx j =

2N 1k =1 k

2N 1j =0 (sin2kx j sin2mx j ) +

4N 1k=2 N k

2N 1j =0 (sin2kx j sin2mx j ) + . . . +

2( p+1) N 1k=2 pN k

2N 1j =0 (sin2kx j sin2mx j ) + . . . .

Insa sin2(2 pN + k)x j = sin 2 kx j pentru x j = j/ (2N ). Se obtine asadarbm = m + m +2 N + m +4 N + . . . ,

altfel spus contribut iile frecvent elor mari m + 2 N, m + 4 N , . . . apar mascate ncoecientul bm . Fenomenul de mascare 8 impune s a urm arim n aplicat ii variat iacoecientilor ak , bk la dublarea lui N . Fie a

(2 N )k , b

(2 N )k coecientii determinati

folosindu-se 2N puncte si a (4 N )k , b(4 N )k coecientii determinati folosindu-se 4N

puncte. Daca a (2 N )k = a(4 N )k , b

(2 N )k = b

(4 N )k pentru k N atunci numarulde puncte a fost considerat sucient de mare pentru a elimina fenomenul de

mascare.Exemplu . Folosind interpolarea trigonometric a s a se aproximeze pozit ia punctelor

situate pe elipsa

x2 /a 2 + y2 /b 2 = 1 .

Rezolvare . Ca n orice problem a de aproximare, alegerea variabilelor este important a.Reprezentarea n coordonate carteziene x, y este dezavantajoas a, conduc and la dou afunct ii

y = b 1 x2 /a 2 .De aceea se prefer a coordonatele polare r, sau coordonatele paramametrice r, t . Incoordonatele polare x = r cos , y = r sin , elipsa este data de funct ia

r () = ab a2 sin 2 + b2 cos 2 1/ 2

,

cu [0, 2]. Funct ia r () are perioada 2 . O aducem la perioada 1 prin transfor-

marea s = / 2 ,

r (s) = ab a2 sin 2 2s + b2 cos 2 2s 1/ 2

.

8 Comportarea este descris a deseori prin denumirea din englez a de aliasing .


29/313


Figura 1.7: Variatia erorii relative la interpolarea trigonometric a a unor elipse. Re-zultatele pentru reprezentarea parametric a sunt unite cu linii.

In reprezentarea parametric a, elipsa este descrisa de x = a cos t, y = bsin t, astfelnc at obt inem

r (t) = a2 sin 2 t + b2 cos 2 t1/ 2

cu t [0, 2]. Aducem funct ia la perioada 1 prin s = t/ 2 si avem o a doua reprezen-tare

r (s) = a2 sin 2 2s + b2 cos 2 2s1/ 2

Vom nota prin gN (s) polinomul trigonometric ce interpoleaz a r (s) n 2 N puncte echi-distant repartizate n intervalul [0 , 1]. Evalu am calitatea interpol arii pentru diversevalori ale lui N = 4 , 8, . . . , 256 prin calculul erorii relative pe o diviziune mai deas a

{j = j/ 2048, j 0, 2048}N =

4N

j =0|r (j ) gN (j )| / |r (j )| .

Presupunem c a a = 1 si vom studia comportarea erorii relative pentru mai multe valoriale lui b, n cele doua reprezentari adoptate .

Rezultatele sunt prezentate n g. 1.7. In toate cazurile, interpolarea este conver-gent a: pe m asur a ce N creste eroarea se aproprie de zero. Cea mai rapid a convergent a

se obtine pentru a/b = 1 / 2 deoarece funct ia ce descrie elipsa are variat ii lente pe inter-valul [0, 1]. Odat a cu scaderea raportului a/b , observ am c a este nevoie de un numarmai mare de termeni n polinomul trigonometric pentru a se obt ine o precizie data.Elipsa are variat ii rapide n punctele ( 1, 0) si este nevoie de mai multi termeni pentrua satisface criteriul Nyquist. In plus, este de asteptat ca aproximarea s a e afectata side fenomenul Gibbs n zona punctelor ( 1, 0) variat iile funct iilor sunt ntr-at at de


30/313


rapide nc at nesatisfacerea criteriului Nyquist le face s a apar a ca niste discontinuit ati.Se poate observa ca, pentru un acelasi raport a/b reprezentarea parametric a (puncteunite cu linii n g. 1.7) conduce la erori mai mici n toate cazurile studiate, arat andimportant a unei alegeri judicioase a reprezent arii.

Sa mai remarcam ca, o dat a atins a o precizie de circa 12 cifre, n cazul unei funct iir (s) netede, a/b = 1 / 2, cresterea lui N nu mai mbun at ateste calitatea aproxim arii.Calculele au fost efectuate n virgul a mobil a cu o precizie de circa 11 , 12 cifre semni-cative. Odata atins acest prag, cresterea lui N produce efecte detrimentale deoareceprin cresterea num arului de operatii aritmetice are loc o acumulare a erorilor de rotu-njire. Dac a ntr-adev ar este nevoie de o acuratet e mai mare trebuie lucrat n preciziesporita.

1.1.4 Interpolare n planul complex

Anumite aplicatii importante impun interpolarea n planul complex. Printreacestea se numar a transformarea conform a si transformata Fourier discret a careau multe aspecte asem anatoare. Transformarea conform a a unui contur ( C ) deforma oarecare din planul complex z = x + iy, pe un cerc (K ) de raz a egalacu unitatea plasat n planul = + i, astfel nc at exteriorul conturului s ase transforme pe exteriorul cercului (g. 1.8), este o problem a des ntalnit an mecanica uidelor incompresibile sau n electricitate. Conform teoremei dereprezentare Riemann, odat a precizate contururile, se mai pot alege trei para-metrii reali. Acestia se determin a de regul a din impunerea unor corespondent ece au semnicat ie zica. Vom ncerca mai nt ai determinarea formei generale atransformarii conforme, l asand la urm a stabilirea celor trei parametrii. Formageneral a a transformarii conforme este o serie care se limiteaza la un num arconvenabil de termeni M 9

z =M 2

n = 1C n n . (1.61)

Cei M coecienti C n , n = 1, 0, . . . , M 2 se determin a din condit ii de cores-pondent a ntre punctele P k (zk ) de pe contur si punctele k ( k ), k = 0 , M 1de pe cerc. Prin urmare se obt ine sistemul de ecuatiizk =

M 2

n = 1C n nk . (1.62)

Este convenabil ca mp art irea cercului s a se faca n part i egale, adic a sa luam

k = exp( ik ), i 1, 2/M . (1.63)Coecient ii C n sunt solut ia sistemului

M 2

n = 1C n e ikn = zk , (1.64)

9 Forma reprezentat a a fost astfel aleas a nc at punctele de la innit ale celor dou a planecomplexe sa corespund a.


31/313


Figura 1.8: Corespondenta punctelor din planul zic si cel de calcul.

punctele zk ind date. Sistemul ( 1.64) se inverseaz a usor, t in and seama c afunct iile exp(

ikn ) formeaz a un sistem ortogonal. Intr-adevar, nmult ind

(1.64) membru cu membru, cu exp( ijk ), j = 1 , 2,...,M 2 si sumand dupa k,se obtine succesivM 1

k=0

zk eijk =M 1

k=0

M 2

n = 1C n eik ( j n ) =

M 2

n = 1C n

M 1

k=0

eik ( j n ) . (1.65)

Dar progresia geometric a de rat ie exp i( j n) se poate suma obtin andu-seM 1

k=0

eik (j n ) =1 ei ( j n )M 1 ei ( j n )

= 0 daca j = nM daca j = n .

In consecint a, din ( 1.65) se obtine

C j =1

M

M 1

k =0

eijk , j 1, M 2 . (1.66)

Trebuie mentionat c a problema nu s-a ncheiat odat a cu gasirea coecient ilorC j pentru un M ales. Anume, mai trebuie vericat a atingerea preciziei doriteprin considerarea a unor puncte intermediare ca, de exemplu,

k = exp [ i(k + 1 / 2)] , k 0, M 1 .Daca axele zk sunt prea ndep artate de conturul ( C ) atunci trebuie ndesit areteaua de puncte de calcul si, n consecint a, marit numarul de termeni n seria(1.61).

Rezultatele obtinute se nscriu ntr-o formulare mai general a, important apentru multe aplicat ii10 . Pentru reteaua de puncte x j = jh , j = 0 , 1, . . . , N 1,10 Algoritmul TFR prezentat n continuare este, conform unor studii statistice de utilizare

a calculatoarelor, al doilea cel mai des utilizat algoritm, ind devansat doar de rezolvareadirecta a sistemelor liniare. Exemple de aplicat ii: modelari meteorologice, analize de vibrat iimecanice, prelucrare de imagini.


32/313


avand pasul h = 1 /N , setul de valori ale unei funct ii u n noduri se noteazau(N ) si se numeste funct ie ret ea . Not and mai departe cu radacina de ordinulN a unit at ii

= exp(2 i/N ), (1.67)

se deneste transformata Fourier discret a direct a (TFDD) a funct iei u(N ) caind o alt a funct ie ret ea, notata U (N ) , ale carei valori n noduri sunt date derelat iile

U (N )k =N 1

j =0u(N )j

jk , k 0, N 1 . (1.68)

Se observa ca (1.66) este un caz particular de TFD. Prescurt am notatia prinadoptarea semnului = pentru TFD dat a de (1.68)

u(N ) = U (N ) .

Sistemul ( 1.68) se poate inversa, prin aplicarea relat iilor de ortogonalitate obt i-nandu-se

u(N )j =1N

N 1

k =0

U (N )k jk , j 0, N 1 . (1.69)

denumita transformata Fourier discret a invers a (TFDI) ce va notata

u(N ) = U (N ) .

Calculul direct al sumelor din ( 1.68) revine la nmultirea vectorului u(N ) cumatricea = [ jk ], j, k 0, N 1 si ar necesita O(N 2) operat ii aritmetice.Insa matricea are propriet at i remarcabile, provenite din structura sa ciclic a,ce permit o evaluare n numai O(N log N ) operat ii. Castigul este enorm, iaralgoritmul ce realizeaza acest c astig se numeste transformata Fourier rapid a sau prescurtat TFR 11 . Vom prezenta cea mai simpl a deducere a algoritmului, oexemplicare a tehnicii generale divide et impera o problema dicila se poaterezolva uneori mai usor prin descompunerea n dou a probleme mai simple.

Vom presupune ca funct ia ret ea u(N ) este denit a pe un num ar par de puncteN = 2 P . Din funct ia ret ea u(2 P ) construim doua noi funct ii ret ea v(P ) , w(P )denite n punctele de indice par, respectiv impar

v(P )m = u(2 P )2m , w

(P )m = u

(2 P )2m +1 , m 0, P

1 .

11 O aplicatie tipica provine din domeniul previziunilor meteo pe termen scurt pe care sebazeaza navigatia aerian a si marin a. Tipic, pentru o previziune pe 3 zile, se calculeaz a 106TFD-uri de lungime N = 2 14 . Un calculator performant, de vitez a 109 operatii aritmetice pesecunda, n-ar reusi calculul pro duselor matrice-vector dec at n 74 ore prea t arziu ca sa maie de folos. Cu folosirea TFR timpul se reduce la circa 4 minute!


33/313


Suma ( 1.68) se poate scrie

U (2 P )k =2P

1

j =0

u(2 P )j jk = P

1

m =0v(P )m

2mk + w(P )m (2 m +1) k =

P 1

m =0v(P )m

2mk + kP 1

m =0w(P )m

2mk . (1.70)

Insa 2 este r ad acina de ordinul N/ 2 = P a unit at ii astfel ncat se observ aaparit ia n ( 1.70) a transformatelor

v(P ) = V (P ) , w(P ) = W (P ) ,

relat iile de legatur a ntre cele trei transformate ind deosebit de simple

U (2 P )k = V (P )

k + k W (P )k , U

(2 P )k+ P = V

(P )k k W

(P )k , k 0, P 1 (1.71)

unde s-a folosit identitatea k + P = k . Acest rezultat este cunoscut ca lema lui Danielson si Lanczos . Evaluarea directa a lui U (2 P ) ar costat O(4P 2)operat ii aritmetice. Acum avem de evaluat dou a transformate, V (P ) , W (P ) ceecare necesit a O(P 2) operat ii pentru un total de O(2P 2) operat ii. Reducereala jum atate a calculelor este benec a, ns a se poate c astiga si mai mult dac a N este o putere a lui 2 , N = 2 q . In acest caz separarea pe indici pari si imparipoate continua p an a cand se ajunge la transformarea unui vector de lungime1 ce nu necesit a nici o operat ie aritmetica deoarece se reduce la transformareade identitate, U (1) = u(1) . Singurele operat ii artimetice ce mai raman suntnmultirile cu k din (1.71). Calculul se poate vizualiza ca parcurgerea arborelui

cu q = log 2N nivele din g. 1.9. Pe ecare nivel sunt necesare doar N operat iiaritmetice astfel nc at obt inem costul mentionat anterior de O(N log2 N ).Separarea repetat a n indici pari si impari conduce la o alt a ordonare a

componentelor vectorului u(N ) . Exemplicam pentru N = 2 3 = 8. Vectorulu(N ) are componente de indici de la 0 la 7. Procedeul de separare par-impar seaplica de dou a ori dup a cum este ar atat n urm atoarea schema

etapa 0: 0 1 2 3 4 5 6 7etapa 1: 0 2 4 6 1 3 5 7etapa 2: 0 4 2 6 1 5 3 7

Inmult irile cu diversele puteri ale lui trebuie s a le ncepem asupra vectoruluipermutat {u0 , u4 , u2 , u6 , u1 , u5 , u3 , u7}. Permutarea se poate construi ns a foarteusor daca scriem indicii init iali si nali n binar

etapa 0: 000 2 0012 0102 0112 1002 1012 1102 1112etapa 2: 000 2 1002 0102 1102 0012 1012 0112 1112

.

Observat i ca inversarea ordinei de citire a indicilor din ultima etap a corespundechiar la numerotarea natural a 0, 1, . . . , 7 anume: 0002 citit de la dreapta la


34/313


Figura 1.9: Arborele de recurent a de la TFR.

st anga este 000 2 = 0, 100 2 citit de la dreapta la st anga este 001 2 = 1, . . . ,0112citit de la dreapta la st anga este 110 2 = 6, 111 2 citit de la dreapta la st angaeste 111 2 = 7. Rezultatul este valabil pentru orice N = 2 q iar permutarea nal aa indicilor se zice n ordine bit-inversat a . Putem da acum algoritmul TFR

j = 1 : N 2k invbit ( j,q)dac a k > j permut a(u j , uk )rad exp( semn 2i/N );0 1 j = 1 : N/ 2

j

rad

j 1

t 1; p 2; s N/ 2

cat timp p N j = 0 : p : N 1r 0k = 0 : t 1tmp u j + ku j + k tmp + r u j + k+ tu j + k+ t tmp r u j + k+ tr r + st 2t; p 2 p; s s/ 2 .

Vectorul initial u este nlocuit de transformata sa discret a ca atare algoritmulnu necesit a memorie suplimentar a. Variabila semn ia valoarea 1 pentru TFRdirect a si 1 pentru TFR invers a; n acest ultim caz vectorul nal mai trebuiemp art it la N . Funct ia invbit ( j,q) ntoarce indicele ce se obtine prin inversareaordinei bit ilor de la 0 la q 1 din j . Operat ia se poate exprima ntr-un limbajde nivel nalt prin nmult iri si mpart iri cu 2, dar de regula este implementat aecient la nivel de cod masina. Vericarea condit iei k > j are rolul de a nustrica permut arile deja efectuate n cadrul buclei.

Algoritmul TFR poate folosit pentru calculul coecient ilor polinomuluitrigonometric ( 1.57) precizat i de relat iile (1.58). O cale evident a este sa se scrie

ck ak + ibk =1N

2N 1

j =0

yj exp(2ijk/ 2N ), k 0, N

si sa se aplice TFR asupra vectorului yj , j 0, 2N 1 ce are 2N componente.Insa, dac a vectorul {yj }este real, se poate obtine rezultatul cu doar jum atate


35/313


Figura 1.10: O masur atoare experimental a a unei vibratii compuse complexe.

din efortul de calcul printr-o separare n componente pare si impare

zm = y2m + iy2m +1 , m 0, N 1 .Introducem transformatele directe z = Z , (y2m ) = P , (y2m +1 ) = I . Avem

Z k = P k + iI k , Nck = P k + I k exp(ik/N ) .

Se verica imediat c a TFD, F a unei funct ii reale f satisface F N k = F k , undeF k este conjugata complex a a lui F k . Ca atare

Z N k = P N k + iI N k = P k + iI k Z N k = P k iI ksi putem deduce expresiile P k = ( Z k + Z N k )/ 2, I k = i(Z k Z N k )/ 2 astfelnc at

ck = 12N (Z k + Z N k ) i2N (Z k Z N k )exp( ik/N ), k 0, N . (1.72)Exist a multe astfel de combinat ii posibile ce furnizeaz a moduri economice de acalcula transformata Fourier a unui vector ce are propriet at i de simetrie supli-mentare.

Exemplu . Prezentam o aplicat ie real a tipic a. Un senzor de vibratii a m asuratdeplas arile unei componente mecanice rezult and datele din g. 1.10. Vibrat iile defrecvent a mai joas a de 1 Hz nu prezentau interes astfel nc at fereastra de timp ncare s-au m asurat datele a fost de t = 1 sec. Se cere identicarea frecvent elor deamplitudine maxim a astfel ncat s a se poat a identica sursele de excitat ie ce producvibrat ia componentei.

Rezolvare . Masur atorile au fost efectuate cu patru rate de esantionare N = 256 ,

512, 1024 rezult and datele u(256)

, u(512)

, u(1024)

. Pentru identicarea frecvent elor do-minante se calculeaza spectrul de putere al semnalului denit de

P u (f ) = 2 U (f )U (f ) ,

unde u = U . Spectrul se calculeaza prin aplicarea TFR asupra datelor din g. 1.10.Deoarece semnalul este real se aplic a relat ia ( 1.72), lungimea transformatelor ind


36/313


Figura 1.11: Spectrele de putere ale semnalului anterior pentru N = 256 , 512, 1024.

N/ 2. Rezultatele sunt prezentate n g. 1.11. Maximele locale sunt denumite picuri si corespund la vibrat iile predominante. Se poate observa c a la trecerea la o rat a deesantionare mai mare uneori apar noi picuri n domeniul frecvent elor mari iar unelepicuri din domeniul frecvent elor mici dispar. De asemenea, amplitudinile asociate unuipic se modic a uneori. Alte picuri nu par afectate de cresterea ratei de esantionare.

Comportarea rezult a din fenomenul de mascare discutat n 1.1.3. De exemplu,picurile , sunt stabile la trecerea de la spectrul 2 la 3. Rezult a ca pentru acestefrecevent e mici criteriul Nyquist este statisf acut si n plus nu are loc mascharea unorfrecvent e mai nalte. Picul din spectrul 2, construit cu N = 512 dispare completns a la trecerea la spectrul 3, construit cu N = 1024. In spectrul 2 picul erafals, n sensul ca nu reprezenta o vibrat ie real a cu frecvent a f . De fapt, rata deesantionare folosit a era prea mica astfel ncat se nregistra contribut ia unei frecventemai nalte f 2f , frecvent a a carei contribut ie devine discernabil a atunci c and sedubleaz a rata de esantionare n spectrul 3. Intr-adevar pentru a discerne frecvent af = 350 Hz criteriul Nyquist indic a necesitatea a cel put in 700 puncte, conditierealizata doar de m asurarea u (1024) . Exemplul arat a important a studiilor de rezolut ien aplicatii de aproximare spectral a, studii n care se urm areste stabilitatea spectruluila modicarea ratelor de esantionare. In aplicat ii practice, fenomenul de mascare seelimin a prin folosirea unor ltre trece-jos cu pret ul pierderii port iunii de frecventenalte a spectrului.


37/313


1.2 Aproximarea mini-max

Trecem acum la considerarea celui de-al doilea criteriu de aproximare ( 1.3).Exemplul lui Runge din 1.1.1 a ar atat ca aproximarea prin interpolare, chiardaca trece prin noduri, poate prezenta erori mari (chiar innit de mari!) ntrenoduri. De aici provine ideea de a ncerca g asirea unei aproximari optimale petot intervalul [ a, b], nu numai n noduri. O exprimare matematic a a acestei ideieste sa se minimizeze eroarea maxima ntre functia f (x) si aproximanta g(x),adic a criteriul ( 1.3) pe care-l repet am aici

maxx [a,b ] |f (x) g(x)| = minim . (1.73)

Astfel formulata problema este momentan incomplet denit a deoarece nu amprecizat din ce multime de funct ii luam pe g(x).

1.2.1 Polinoamele Cebasev

Vomncepe prin a presupune iar asi c a funct ia aproximant a g(x) apart ine mult imiipolinoamelor. Criteriul de minimizare a erorii maxime ( 1.73) face referire laun anume interval [ a, b]. Deoarece rezultatele de mai jos depind de interval,adopt am un interval canonic [ 1, 1]. O transformare de scara x = z(ba)/ 2 +(b+ a)/ 2, x [a, b], z [1, 1] poate aduce orice alt interval nit pe cel canonic.De data aceasta nu vom scrie funct ia aproximant a g(x) ca o combinat ie li-niar a de monoame. Examinnd comportarea monoamelor xk pe intervalul [ 1, 1]se constata ca toate iau valori absolute maxime la capetele intervalului; prin ur-mare este de asteptat ca erorile de aproximare s a nu e distribuite uniform peinterval. Aceasta sugereaz a cautarea unor polinoame care s a nu varieze mono-ton pe [

1, 1]. Intuitiv, variatia monotona face dicila descrierea unor variat ii

rapide ale funct iei f n interiorul intervalului [ 1, 1]. Ceea ce ne trebuie este omult ime de polinoame denite pe [ 1, 1] care sa poat a descrie astfel de variatiintr-un mod c at mai economic, adic a folosind un num ar mic de polinoame.Funct iile trigonometrice considerate anterior prezentau part ial o astfel de com-portare: cos2 x si sin 2x descriau variatiile de perioad a 1, cos4x si sin 4xdescriau variatiile de perioad a 1/ 2, s.a.m.d. Spunem part ial deoarece nu a fostndeplinit a si dorint a de economicitate avem dou a funct ii, sin si cos pentrudescrierea variatiilor de o anumita perioad a. Pe baza funct iilor trigonometriceputem ns a introduce o clasa remarcabila de polinoame descoperite de Ceb aseva caror proprietat i sunt deosebit de favorabile aproxim arii optimale c autate naceast a sectiune.

Vom lua z [1, 1]. Introducem variabila [0, ] prin relat iilez = cos , = arccos z .

Sa consider am acum funct ia T n (z) = cos n . Funct ia are n rad acini pe (1, 1)T n (z) = cos n = 0 k = (2 k 1)/ 2n, k 1, n ,


38/313

1.2. Aproximarea mini-max 31

Figura 1.12: Primele 5 polinoame Ceb asev.

Tabelul 1.2: Polinoame Cebasevn T n (z) zn

0 1 T 01 z T 12 1 + 2z2 (T 0 + T 2)/ 23 3z + 4 z3 (3T 1 + T 3)/ 44 18z2 + 8 z4 (3T 0 + 4 T 2 + T 4)/ 85 5z 20z3 + 16 z5 10T 1 + 5 T 3 + T 5

zk = cos k = cos[(2 k 1)/ 2n], k 1, n .Din identitatea trigonometric a cos n +cos( n2) = 2 cos( n1) cos deducemo relatie de recurent a pentru T n (z)

T n (z) = 2 zT n 1(z) T n 2(z) . (1.74)Deoarece T 0(z) = 1 , T 1(z) = z iar T 2(z), T 3(z), . . . rezult a din ( 1.74) prinoperat ii de adunare si nmult ire, putem observa ca T n (z) sunt polinoame degradul n, denumite polinoame Ceb asev relative la intervalul compact [ 1, 1].Funct ia cos n se numeste funct ia generatoare a acestor polinoame.

Gracele primelor c ateva polinoame Cebasev sunt prezentate n g. 1.12 iarexpresiile lor n tabelul 1.2 mpreun a cu exprim arile monoamelor zn n funct iede T

k(z), k 0, n . Se observa din g. 1.12 ca ecare polinom Cebasev de grad

mai nalt poate descrie variat ii mai rapide ale unei functii pe intervalul [ 1, 1].

Din formula de recurenta, se observa ca n T n (z) coecientul lui zn este2n 1 , astfel nc at polinomul T n (z) 21 n T n (z) are coecientul lui zn egal


39/313


cu unitatea. Polinoamele cu un coecient 1 al termenului de grad maxim senumesc polinoame monice . Proprietatea esent ial a a polinoamelor Cebasev ce leface adecvate aproximarii mini-max este:

Teorema. Dintre toate polinoamele monice de grad n xat, T n (z) are cea maimica margine n valoare absolut a pe intervalul [ 1, 1].

Demonstrat ie. Proced am prin reducere la absurd. In acest scop s a observ ammai nt ai ca funct ia T n (z) = 2 1 n |cos n| ia de n + 1 ori valorea maxima 21 npe [1, 1], anumen punctele distincte zk = cos( k/n ), k 0, n . Sa consider amapoi un alt polinom monic de acelasi grad pn (z), si sa presupunem prin absurdca acesta ar avea o margine superioar a mai mica n modul dect T n (z) pe [1, 1],adic a

supz [

1,1]

| pn (z)| < supz [

1,1]

T n (z) = 2 1 n . (1.75)

In consecint a, chiar n punctele zk

| pn (zk )| < 21 n . (1.76)Considernd diferent a dn 1(z) a polinoamelor T n (z) si pn (z), dn 1(z) T n (z) pn (z) vom obt ine, evident, un polinom de gradul n 1. Din (1.76) rezult a

(1)k T n (zk ) pn (zk ) = ( 1)k + n dn 1(zk ) > 0, k n, 0adica dn 1(z) are n schimb ari de semn pe (1, 1), deci polinomul de graduln 1, dn 1(z), ar avea n rad acini, ceea ce este absurd. Rezult a ca presupunerea(1.75) nu este adev arat a, deci oricare ar polinomul monic pn (z) avem

supz [ 1,1] | pn (z)| 2

1 n . 2 (1.77)

In continuare se dau dou a aplicat ii importante ale polinoamelor mini-max.

1.2.2 Minimizarea erorii la interpolarea polinomiala

Trecem la un prim exemplu concret de aplicare a criteriului mini-max, anumenformularea generala (1.73) vom lua funct ia aproximant a un polinom de gradulN 1

maxx [a,b ] |f (x) pN 1(x)| = minim .

Dorim s a determinam, dintre toate polinoamele de grad N 1, pe cel ce mini-mizeaza abaterea maxim a fat a de funct ia f (x). Polinomul obt inut va numitpolinom mini-max de grad N 1 al functiei f (x) pe intervalul [ a, b]. Este conve-nabil s a privim polinomul pN 1 ca ind denit de faptul ca trece prin nodurile{(x i , yi ), i 1, N }. Nodurile nsa, spre deosebire de problema anterioar a de


40/313

1.2. Aproximarea mini-max 33

interpolare, sunt acuma necunoscute . Vom minimiza eroarea maxima printr-oalegere adecvata a nodurilor. Eroarea este n acest caz chiar restul interpol ariice are forma ( 1.16)

RN 1(x) =N

i =1

(x x i ) f (N ) ()/N ! .

In general derivata f (N ) () nu este cunoscuta, astfel ncat se poate pune doarproblema minimizarii produsului. Acest este un polinom de gradul N . Utilizandschimbarea de variabil a

x = z(ba)/ 2 + ( b + a)/ 2 (1.78)trecem de la intervalul [ a, b] pe care dorim sa minimiz am eroarea la intervalul[1, 1]. Se obtine, consider and si ( 1.77)

N

i =1(x x i ) = ba2

N N

i =1(z zi ) ba2

N

21 N .

Rezult a ca optimul, corespunzand egalit at ii, se realizeaz a dac a punctele zi suntrad acinile polinomului Cebasev de gradul n. Vom aranja indicii astfel ncatrad acinile sa apar a n ordine crescatoare

zN i +1 = cos[(2 i 1)/ 2N ], i 1, N . (1.79)In aceste condit ii se obt ine cea mai mic a margine superioara pentru restul

la interpolare

|RN 1(x)| 21 2N (ba)N max [a,b ] f (N ) () /N ! .

Marginea depinde de marimea intervalului ( a, b), de numarul de noduri N si de derivata f (N ) . Rezultatul obtinut se poate enunt a: dintre toate poli-noamele de interpolare de grad N 1, cel ce minimizeaza eroarea maximamax |f (x) pN 1(x)| cu x [1, 1] este cel construit cu abscisele nodurilordate de r ad acinile polinomului Cebasev de grad N .

Exemplu . Am vazut c a interpolarea funct iei lui Runge f (x) = 1 / (1 + x2 ) pe[5, 5] cu o repartit ie echidistant a a nodurilor nu converge. Vom determina acum unalt polinom ce aproximeaz a pe f (x) dar la care nodurile nu mai sunt luate echidistant,ci sunt determinate de criteriul mini-max, anume vom lua xk = 5 zk , k 1, N , cu zkdeterminate de ( 1.79). Ca n exemplele precedente calitatea aproxim arii este apreciat aprin evaluarea erorii relative = 500i =1 |f (u i ) pN 1 (u i )| / |f (u i )| pe o diviziune maina a intervalului [ 5, 5]. Variatia erorii cu N este

N 5 10 15 20 25 30 35 40 45 50lg

0.4

0.8

1.2

1.6

2.1

2.4

3.0

3.3

3.8

4.2

oberv andu-se c a aproximarea converge. Observ am acum clar important a alegerii no-durilor pe care se bazeaza o interpolare. Aproximarea prezent a este tot o interpolare,bazata ns a pe alte noduri decat cele echidistante considerate anterior. Comportareaaproxim arii este redata n g. 1.13. Se observ a ca, pentru N mic, apar oscilat ii ntrenoduri dar, spre deosebire de cazul redat n g. 1.1, acestea nu mai cresc la innit.


41/313


Figura 1.13: Polinoamele de interpolare cu noduri alese optimal pentru N = 11 , 21, 31suprapuse peste f (x). Nodurile sunt reprezentate doar pentru N = 31.

1.2.3 Aproximarea aproape mini-max a unei funct ii

Vom considera acum o relaxare a criteriului ( 1.73) anume

maxx [a,b ] |f (x) pN (x)| e . (1.80)

Criteriul ( 1.80) este mai larg dec at cel anterior ( 1.73) deoarece nu mai impunemdeterminarea polinomului aproximant pN (x) astfel nc at s a realiz am un minimal erorii, ci doar ca eroarea s a devin a mai mic a decat un prag dat e. De dataaceasta gradul polinomului nu se mai consider a xat, ci c aut am polinomul degradul cel mai mic ce satisface conditia ( 1.80). Asemenea probleme apar tipic nsituat ii n care dorim o aproximare c at mai economic a n operatii aritmetice 12 .

Vom lua [1, 1] pentru intervalul din ( 1.80) si vom presupune ca funct ia f (x)are o dezvoltare n serie de puteri

f (x) =M

k =0

bk xk , (1.81)

unde eventual putem avea M . Puterile xk se pot exprima n funct ie depolinoamele Ceb asev (vezi tabel 1.2) astfel nc at f se poate scrie

f (x) =

M

k=0ak T k (x) . (1.82)

Am vazut c a polinoamele Ceb asev sunt mai eciente n descrierea variat iilorunei funct ii pe [1, 1] decat monoamele xk . Este asadar de asteptat ca sirul

12 Procedeul ce urmeaz a mai este cunoscut sub denumirea de economizare de serii .


42/313


43/313


Figura 1.14: Variatia erorii e(z) = cos z p6 (z).

ce realizeaz a eroarea impusa dar cu 40% mai putine operatii aritmetice dec at S 10 .Gracul erorii cos x p6 (x) prezentat n gura ( 1.14) conrm a atingerea precizieiimpuse.

1.3 Aproximarean sensul celor mai mici p atrate

Reamintim ca, n acest caz, criteriul de aproximare l reprezint a minimizareasumei (1.5)

S =N

k=1

[yk g(xk )]2 = minim, (1.83)

valorile (xk , yk ) ind date. Este convenabil ca aproximanta g(x) sa se pun a subforma (1.6)

g(x) =n

j =1

a j gj (x) , (1.84)

gj (x) ind funct ii cunoscute, liniar independente, iar a j , j = 1 , 2, . . . , n parame-tri nedeterminat i. Ca si n 1.2, criteriul ( 1.83) pune o problem a de minimizarea erorii. Vom vedea cum c ateva not iuni mai avansate, considerate n 1.4, ne

vor conduce la concluzia c a at at criteriul mini-max c at si cel al celor mai micipatrate sunt doar formul ari diferite ale unui acelasi criteriu de aproximare op-timal a. Deocamdata vom considera nsa doar aspectele simple ale teoriei.

Aproximarean sensul celor mai mici p atrate este utilizat a mai ales n cazulprelucr arii datelor experimentale. In acest caz, nu se recomanda folosirea inter-polarii deoarece valorile m asurate contin erori inerente, repartizate probabilistic


44/313

1.3. Aproximarea n sensul celor mai mici p atrate 37

avand caracterul unor perturbat ii care trebuie, dimpotriv a, eliminate 13 . Estedeci evident c a num arul de parametri n trebuie s a e mai mic decat numarulde noduri N, n < N .

Criteriul ( 1.83), ca si 1.73, minimizeaz a o eroare. In acest caz nsa eroareaeste scris a ca o suma de p atrate. Consecint a esent ial a a acestei alegeri esteca S (a1 , . . . , a n ) este derivabil a, fapt ce poate exploatat pentru a determinacoecientii a i . Pentru ca S (a1 , . . . , a n ) sa aibe un extrem trebuie ca derivatelepart iale n raport cu a i sa se anuleze

S/a i = 0 , i 1, n. (1.85)

Observati introducerea unui indice i diferit de cel de nsumare j din (1.84).Relat iile (1.85) reprezinta un sistem de n ecuat ii cu n necunoscute. DeoareceS este o sum a de p atrate extremul dat de condit iile (1.86) exist a si reprezintachiar un minim. Folosind expresiile ( 1.83) si aranjand termenii se obtine

n

j =1

a jN

k =1

gi (xk )gj (xk ) =N

k=1

yk gi (xk ), i 1, n , (1.86)

adica un sistem de n ecuat ii pentru cei n parametrii a j , j 1, n .Forma sumei S se poate generaliza ntruc atva introduc and o funct ie pondere

w(x), pozitiv a, continua, cu valori n intervalul [0 , 1], care sa ia n considerareunele distinct ii privind important a valorilor luate n noduri. Se va scrie atunci

S =N

k =1

w(xk )[yk g(xk )]2 = minim, (1.87)

iar sistemul ( 1.86) se va nlocui cun

j =1

a jN

k=1

w(xk )gi (xk )gj (xk ) =N

k=1

yk w(xk )gi (xk ), i 1, n . (1.88)

Sistemele ( 1.86) sau (1.88) se pot rezolva prin metodele din capitolul 5, derezolvare a sistemelor de ecuatii liniare. O observatie important a este legat a defaptul c a aceste sisteme pot ridica probleme legate de introducerea unor erorimari n calculul numeric, mai ales c and diviziunile sunt egale 14 . Funct iile gj (x)se pot alege din baza canonica

gj (x) = x j 1 , j 1, n , (1.89)

sau baze formate din diverse polinoame ortogonale pe mult imea discreta depuncte xk , k 1, N (vezi 1.4.2). Alegerea de polinoame ortogonale are avantajulconsiderabil al reducerii erorilor ce pot apare la rezolvarea sistemelor ( 1.86) sau

13 Un exemplu ar eliminarea zgomotului de fond de la aparatura radio.14 In limbajul din capitolul 5, matricea sistemului ( 1.86) sau ( 1.88) este r au condit ionat a ,

iar rezolvarea cere tehnici speciale cum ar descompunerea n valori singulare.


45/313


Figura 1.15: Abaterile geometrice e (n microni) nregistrate la prelucrarea mecanic aa 1000 de piese.

(1.88). Deseori, din informat ii suplimentare, cunoastem forma cea mai adecvat aa funct iilor gj (x). Un exemplu este prezentat la sf arsitul sectiunii. Folosireaunei combinat ii liniare ( 1.84) pentru functia aproximant a g(x) este convenabiladeoarece conduce la un sistem liniar de ecuat ii pentru coecientii a i . Se potadopta ns a si alte forme, cu complicatia posibilei aparitii a unui sistem neliniarde ecuat ii, mult mai dicil de rezolvat. Incheiem cu observatia c a nu s-a precizatnc a ce grad al polinomului ar conduce la cea mai mic a eroare, adic a la cea maimica valoare a lui S n ( 1.83). Un criteriu util pentru alegerea gradului esteminimizarea expresiei

U = S/ (N n) (1.90)prin varierea lui n .

Exemplul 1 . In cadrul unei operat ii de prelucrare mecanic a se nregistreaz a aba-terile din g. 1.15 de la cotele nominale. Sunt reprezentate si marginile ce determin arebuturi. Se pune ntrebarea dac a are loc vreo crestere sistematic a a abaterilor datorit auzurii utilajului.

Rezolvare . Construim o aproximare liniar a prin cele mai mici p atrate. Procedeuleste denumit regresie liniar a . Expresia pentru abaterile e funct ie de numarul de pieseprelucrate p este

e = ap + b

iar coecient ii a, b rezult a din rezolvarea sistemului

N k =1 p

2k a +

N k =1 pk b =

N k =1 ek pk

N k =1 pk a + Nb = N k =1 ek

Se obtine a = 1 .11 10 2 / piesa si b = 4.78 ceea ce indica o crestere cu 1 aabaterii la ecare 100 de piese prelucrate, crestere datorat a probabil uzurii utilajului.De asemenea din faptul c a reglajul init ial a fost efectuat astfel nc at s a se produc a oabatere negativ a de circa 5 se poate deduce ca aparit ia uzurii era prev azut a.


46/313

1.3. Aproximarea n sensul celor mai mici p atrate 39

Exemplul 2 . O serie de masur atori chimice au furnizat ratele de react ie din tabelulde mai jos. Din teoria react iilor chimice se cunoaste dependent a ratei de react ie de

temperatur a

k = cT exp(K/ RT ) .Se cere estimarea parametrilor c,,K pe baza m asur atorilor efectuate. Constantagazelor R=8314 J/mol/K este cunoscut a.

T k T k T kK mol/cm 3 K mol/cm 3 K mol/cm 3

2000 8.4765e12 2600 5.0189e12 3200 3.3146e122100 7.6895e12 2700 4.6544e12 3300 3.1169e122200 7.0071e12 2800 4.3282e12 3400 2.9364e122300 6.4118e12 2900 4.0351e12 3500 2.7711e12

2400 5.8892e12 3000 3.7708e12 3600 2.6194e122500 5.4280e12 3100 3.5317e12 3700 2.4799e12

Rezolvare . Dependenta k(c,,K ) este neliniara. Putem ns a logaritma relat ia demai sus

ln k = ln c + ln T K/ RT si reobt inem o dependent a liniar a de parametrii c ,,K . Se construieste suma

S =N

j =1

(ln c + ln T j K/ RT j ln kj )2 .

Condit iile de extrem S/ (ln c) = 0, S/ = 0, S/K = 0 conduc la sistemul

N ln c + N j =1 ln T j N j =1 1/T j K/ R = N j =1 ln kjN j =1 ln T j ln c +

N j =1 ln

2 T j N j =1 ln T j /T j K/ R = N j =1 ln kj ln T jN j =1 1/T j ln c +

N j =1 ln T j /T j N j =1 1/T 2j K/ R = N j =1 ln kj /T j

a carui solut ie este ln c = 44 .97, c = 3 .39 1019 mol/K 2 /cm 3 /s, = 2, K = 46000J/mol .Exemplul 3 . Sa presupunem date m asur atori zice ale energiei de vibrat ie ale

unei molecule de O 2 . Se cunoaste dependent a energiei de temperatur a

e = hv/kT exp( hv/kT ) 1

RT

Se cere , frecvent a fotonilor emisi la saltul ntre dou a nivele energetice.Rezolvare . In acest caz nu avem nici o posibilitate de a liniariza problema. Not am

x = hv/k . Suma patratelor este

S =N

j =1

x/T jexp( x/T j ) 1

RT j ej2

,


47/313


Tabelul 1.3: Coecient ii polinomului aproximant din exemplul 4.

a1 a2 a3 a4 a5 a6 a7Simpla -526.0 4769 4059 1242 -1683 11.79 0.7434Dubl a 40.00 10.00 5.000 3.000 2.000 1.000 1.000

iar din S/v = 0 se obt ine o ecuat ie neliniara de forma f (x) = 0 , sau explicit

N

j =1

x/T jexp( x/T j ) 1

RT j ej[exp(x/T j ) 1] x/T j exp( x/T j )

[exp(x/T j ) 1]2= 0

Determinarea solut iei se poate efectua prin metodele din capitolul 4.Exemplul 4 . (Dorn) Dam un exemplu al erorilor mari ce pot ap area atunci cand

nu lucr am cu funct ii ortogonale. Gener am un set de noduri

{xk , yk , k 1, 15

}unde

luam xk = k 1, yk = f (xk ) cu f (x) = 40 + 10 x + 5 x2 + 3 x3 + 2 x4 + x5 + x6 . Vomncerca reobt inerea coecientilor polinomului f prin metoda celor mai mici p atratefolosind alegerea ( 1.89), gj = x j 1 , j 1, 7. Coecient ii obt inut i n urma rezolv ariisistemului liniar 15 n simpla (32 biti) si dubla precizie (64 bit i) sunt prezentat i n tabel1.3. Acumularea rezultatelor part iale s-a f acut n precizie extins a (80 biti) n ambelecazuri. Se poate observa c a obtinem abateri mari de la valorile exacte ale coecient ilorn simpla precizie. M arirea preciziei furnizeaz a rezultatul exact, astfel nc at suntemasigurat i ca modul de calcul este corect si abaterile din calculul n simpla precizie t inde ns asi natura problemei. Exemplul este reluat la sf arsitul sectiunii 1.4 cu folosireaunor funct ii ortogonale pentru gj pentru a vedea efectul benec al acestora.

Sa presupunem acum c a nu am cunoaste gradul polinomului ce ar conduce la ocea mai mic a eroare. Rezolv and problema pentru diverse valori ale lui n se obtinurm atoarele valori pentru S/ (N

n)

n 1 2 3 4 5 6 7 8 9S/ (N n) 5.3E 11 5.1E 10 1.8E 09 1.3E 07 9.4E 05

Met Ode Nu Me Rice

Documents

Transcript of Met Ode Nu Me Rice