Criterii de Aproximare a Functiilor. Aproximarea Functiilor In Sensul Lui Cebisev

1. CRITERII DE APROXIMARE A FUNCȚIILOR

1.1. Introducere

În foarte multe aplicații practice apare necesitatea aproximării unei funcții f:a,b→R printr-o alta funcție F(x) relativ simplă, astfel ca pentru orice valoare a lui x, valoarea lui F(x) sa fie "suficient de aproape" de valoarea lui f(x).

Vom scrie f(x)≈F(x), x € a, b.

Există în special două cazuri în care se impune aproximarea funcției f(x). Primul este acela în care funcția f(x) are o expresie complicata sau este dificila de evaluat sau de manipulat în calcule. Astfel, de exemplu, pentru evaluarea funcției cos(x) prin operații aritmetice se impune mai întâi aproximarea funcției printr-o suma parțiala a seriei de puteri: cos x ≈ 1 – x2 + x4 – … (-1)n x2n

2! 4! (2n)!

Al doilea caz în care se impune aproximarea funcției f(x) este acela în care aceasta este data printr-o tabela de valori, obținuta, de exemplu, ca urmare a unor măsurători:

În această situație se aproximează funcția data tabelar printr-o expresie analitica, care sa permită interpolarea în tabela de valori, cu alte cuvinte estimarea valorilor f(x) pentru x ≠ xk.

Fie M – { f / f: a,b → R } un spațiu vectorial si fie o mulțime de funcții φ0(x), φ1(x), …φk(x), aparținând lui M, liniar independente, adică c0φ0(x)+ c1φ1(x)+….+ ckφk(x)=0, sa rezulte c0= c1=………=ck= 0.

Aproximarea unei funcții f oarecare din M se face printr-o combinație liniara de un număr finit m de funcții de tipul φk, adica f(x) ≈ Fm(x), unde Fm(x) = c0φ0(x)+ c1φ1(x)+….+ cmφm(x) = ckφk(x). Vom numi funcția Fm(x) polinom generalizat, aproximarea funcției f făcându-se în acest caz prin polinoame generalizate.

Foarte frecvent în procesul de aproximare se iau drept set de funcții liniar independente, funcțiile 1, x, x2, …..xm. În acest caz, polinomul de aproximare Fm (x) va fi un polinom algebric. Polinoamele sunt ușor de evaluat, iar suma, diferența si produsul a doua polinoame conduc de asemenea la polinoame. În plus, polinoamele pot fi derivate si integrate cu ușurința. Aproximarea polinomiala se bazează pe teorema de aproximare a lui Weierstrass care arata ca daca f(x) este continua pe intervalul închis a,b atunci pentru orice ε > 0, exista un polinom pn(x) de gradul n=n(ε) , astfel ca: f(x)-pnx< ε, a ≤ x ≤ b.

Din nefericire criteriile existente pentru generarea polinomului de aproximare nu garantează în nici un fel că polinomul găsit este cel pus în evidență de teorema lui Weierstrass.

Un alt set de funcții liniar independente, des utilizate în teoria aproximării, sunt: ½, cos x, sin x cos 2x, sin 2x, … cos mx, sin mx. În acest caz polinomul de aproximare poartă numele de polinom trigonometric.

Fm(x) = a0/2 + a1 cos x + b1 sinx + …+ am cos mx + bm sin mx = a0/2 + (ak cos kx + bk sin kx).

Din expresia funcției Fm(x) se observa ca nu este suficienta cunoașterea funcțiilor liniar independente φk(x), fiind necesara de asemenea determinarea coeficienților ck.. Pentru calculul acestor coeficienți sa presupunem ca spațiul M se poate organiza ca un spațiu metric, adica putem defini pe M o funcție ce măsoară distanta dintre doua funcții oarecare f si g. Vom determina polinomul generalizat Fm(x) deci si coeficienții c0, c1,………,cm impunând condiția ca distanta dintre funcția data f si mulțimea polinoamelor generalizate sa fie cit mai mica. În funcție de modul de definire a distantei se pot pune în evidenta următoarele criterii mai des utilizate

1.2. Criteriul de aproximare prin interpolare

Daca în spațiul metric M se definește distanta dintre doua puncte f si g în felul următor; d(f,g) = f(xi) – g(xi), se ajunge la criteriul de aproximare prin interpolare. Daca se înlocuiește funcția g(x) prin polinomul de interpolare Fm(x), condiția ca distanta dintre funcția f(x) si polinomul de interpolare Fm(x) sa fie minima impune sistemul de ecuații: Fm(xi) = f(xi), i=0, 1, 2…., n.

Daca functia f este data printr-o tabla de valori:

Xk x0 x1……xn

f(xk) f(x0) f(x1)….f(xn)

alegând valorile xi ca fiind cele din tabelă, sistemul de ecuații impune de fapt coincidenta dintre funcție si polinom în punctele xi numite noduri de interpolare (fig.1).

Fig. 1

Se observa din figura ca intre doua noduri de interpolare abaterea polinomului de interpolare de la funcția data ∆(x) = f(x) – Fm(x) poate atinge valori apreciabile. Coeficienții ck ai polinomului de interpolare dat astfel:

Fn(x) = c0 +c1x + c2x2 +… + cnxn

Se calculează conform celor arătate mai sus, rezolvând sistemul de n+1 ecuații liniare cu (n+1) necunoscute c0, c1, …, cn.

c0 +c1x0 + c2x02 +… + cnxn0 = f(x0)

c0 +c1x1 + c2x12 +… + cnxn1 = f(x1)

– – – – – – – – – – – – – – – – – – – – – – – –

c0 +c1xn + c2xn2 +… + cnxnn = f(xn)

Determinantul sistemului

este de fapt determinantul Vandermonde si acesta este diferit de zero daca xi≠xj, i≠j, deci dacă nodurile de interpolare sunt distincte. Rezultă un sistem (c0, c1, ……cn) unic, deci polinomul de interpolare este unic. Schimbând poziția nodurilor de interpolare se poate micșora valoarea abaterii maxime de la funcția dată.

Este recomandabil să se ia nodurile de interpolare dispuse analog zerourilor polinomului Cebîșev de același grad cu polinomul de aproximare.

Deci pentru polinomul de interpolare de gradul n nodurile de interpolare vor fi calculate cu formula:

Se poate pune problema daca mărind gradul polinomului de interpolare (deci făcând pe n sa tindă spre infinit) șirul corespunzător de polinoame pn(x) converge către f(x).

Pentru funcții f(x) relativ simple și pentru anumite alegeri ale nodurilor de interpolare (de exemplu ca fiind zerourile polinoamelor Cebîșev) convergența poate avea loc.

În general lucrând cu noduri de interpolare egal distanțate funcția eroare ∆n(x) = f(x) – pn(x) are ca reprezentare grafică o curbă oscilantă, așa cum apare în fig. 2.

Fig. 2

În interiorul intervalului de convergenta funcția eroare tinde spre zero atunci când n tinde spre infinit, dar în exteriorul acestui interval poate să crească oricât de mult. În afara interpolării descrise mai sus, cu ansambluri de noduri simple, se folosește interpolarea cu noduri multiple, în special duble. Într-un nod multiplu de interpolare, de ordinul k+1, se impune nu numai coincidența valorilor funcțiilor f(x) cu Fn(x), ci și a valorilor derivatelor lor până la ordinul k. Deci, în cazul interpolării cu ansambluri de noduri de interpolare multiple, sistemul de ecuații ce trebuie satisfăcut de către funcția de abatere este:

∆n(xk) = f(xk) – F n (xk) = 0

∆'n(xk) = 0

∆'n(xk) = 0

–––––

∆ n (m)(xk) = 0

Fiecare nod multiplu de interpolare de ordinul (k+1) se poate considera ca o coincidență a k+1 centrii simpli în cazul tendinței spre zero a distantei între aceste centre. De aceea fiecare nod multiplu de ordinul k+1 trebuie socotit drept k+1 noduri simple. Geometric, în cazul unui nod dublu, curba funcției abatere este tangentă la axa Ox în abscisa nodului luat drept considerare. De remarcat că există o legătură între numărul punctelor în care funcția data este tabelata si gradul polinomului de interpolare. Din aceasta cauza o tabela extinsa conduce la un polinom de aproximare de grad mare, ceea ce constituie un dezavantaj. Interpolarea, ce folosește tot intervalul dintre nodurile extreme ale tabelei de valori, poarta numele de interpolare globala.

S-a observat ca mărind gradul polinomului de interpolare globală, eroarea nu se micșorează întotdeauna, ba chiar adesea crește foarte mult pe unele porțiuni ale intervalului de interpolare. O soluție pentru remedierea acestei situații consta în renunțarea la interpolarea globală a funcției și trecerea la interpolarea pe porțiuni sau, altfel spus, la utilizarea funcțiilor spline. O funcție spline polinomiala de gradul k (k≥1) relativ la diviziunea cu n noduri interioare: a=x0 < x1 < x2<…… <xn+1 = b este prin definiție o funcție S(x) ce îndeplinește următoarele condiții:

1) pe fiecare subinterval [ xi, xi+1], 0 ≤ i ≤ n, S(x) coincide cu un polinom de gradul k cu coeficienți reali;

2) S(x) este o funcție de clasa Ck-1[a,b] . Se poate verifica că mulțimea funcțiilor spline de gradul k este un spațiu vectorial de dimensiune n + k + 1, având ca bază funcțiile:

unde:

Rezultă că pentru determinarea unei funcții spline de gradul k se impun n + k + 1 condiții.

În general dacă f : a,b R este o funcție de clasa , atunci exista o funcție spline S(x) unica, de grad 2k-1 ce aproximează pe f(x) si astfel ca următoarele n+2k condiții să fie verificate

S(x1) = f(x1), ……..S(xn) = f(xn)

S(a) = f(a), S'(a) = f'(a), …… S(k-1)(a) = f(k-1)(a)

S(b) = f(b), S(b) = f'(b), …… S(k-1)(b) = f(k-1)(b)

Aplicația 1

Sa se determine funcția spline de grad 1 pe intervalul [0, π] asociata funcției f(x) = cos x si nodurilor x1 = π/3, x2 = π/2. În acest caz 1, x, (x- π/3)+, (x- π/2)+

constituie o baza si funcția S(x) devine:

S(x) = C1 + C2x + C3(x- π/3) + + C4 (x- π/2)+

Pentru determinarea celor 4 coeficienți, condițiile de mai sus devin:

S(π/3) = cos π/3 = 1/2

S(π/2) = cos π/2 = 0

S(0) = cos 0 = 1

S(π) = cos π = -1

Obtinem sistemul:

C1 + C2π/3 = 1/2

C1 + C2π/2 + C3π/6 = 0

C1 = 1

C1 + C2π + C3 2π/3+ C4 π/2 = -1

Rezulta:

C1 = 1, C2 = -3/2π, C3 = -3/2π, C4 = 1/π

si deci

cos x ≈ S(x) =

Daca în general prin interpolare se intelege estimarea valorilor unei functii y(x) date tabelar pentru valori ale argumentului x situat intre nodurile de interpolare, atunci interpolarea inversa semnifica, asa cum este de asteptat, aproximarea argumentului x ce corespunde la o valoare y(x) data, presupunand cunoscut sirul valorilor y, y1, ….yn. Este evident ca în interpolarea inversa se schimba rolul variabilelor xk si yk. În plus suntem în situația în care noile argumente yk nu sunt egal distanțate. Prin extrapolare se înțelege estimarea valorilor funcției y(x) în puncte care se găsesc în afara intervalului în care cad argumentele date x0, x1, ……xn

1.3. Criteriul de aproximare cu abatere medie pătratica minima

Presupunem ca suntem în spațiul L2 [a,b] al funcțiilor de pătrat integrabil pe [a,b]. Dacă distanta dintre funcțiile f si g este data prin:

Atunci aproximarea ce folosește această distanță poartă numele de aproximare cu abatere medie pătratică ponderată minimă. Atunci α(x) > 0 este o funcție numită pondere (în particular α(x) = 1).

Adeseori se definește distanța

unde x0, x1, ……xn sunt n+1 puncte cunoscute din intervalul (a,b). Dacă în expresia distanței de mai sus funcția g(x) se înlocuiește cu Fm(x), definită ca un polinom algebric de grad m ≤ n, se ajunge tot la o aproximare cu abatere medie pătratica minimă, metoda de aproximare purtând denumirea de metoda celor mai mici pătrate a lui Gauss. Această metodă este frecvent utilizată în prelucrarea matematică a datelor experimentale. Metoda celor mai mici pătrate se utilizează ori de câte ori exista îndoieli privind exactitatea valorilor f(x) culese experimental. În acest caz nu este recomandabil sa se impună identificarea polinomului de aproximare cu funcția data în valorile xi așa cum se întâmplă în cazul interpolării. Deseori se întâmplă ca polinomul de aproximare dorit sa aibă grad mic, dar exista la dispoziție un mare număr de date xi. Deoarece criteriul de interpolare face o legătura intre gradul polinomului de aproximare si numărul de valori xi., se recomanda în acest caz metoda celor mai mici pătrate care nu impune o asemenea condiție. După cum s-a văzut, aproximarea medie pătratica cere ca polinomul de aproximare sa se calculeze astfel ca valoarea abaterii

sa fie minima (fig. 3)

Fig. 3

Din punct de vedere geometric, metoda aproximării cu abatere medie pătratica minima conduce la cerința ca aria porțiunii hașurate sa fie minima.

1.4. Criteriul de aproximare în sensul lui Cebîsev

Presupunem în acest caz ca spațiul metric M este C[a,b], adică mulțimea funcțiilor continue definite pe [a,b] distanta dintre doua funcții f si g fiind data de :

Daca în locul funcție g(x) punem polinomul de aproximare Fm(x), m≤ n calculul coeficienților ck se va face din cdentificarea polinomului de aproximare cu funcția data în valorile xi așa cum se întâmplă în cazul interpolării. Deseori se întâmplă ca polinomul de aproximare dorit sa aibă grad mic, dar exista la dispoziție un mare număr de date xi. Deoarece criteriul de interpolare face o legătura intre gradul polinomului de aproximare si numărul de valori xi., se recomanda în acest caz metoda celor mai mici pătrate care nu impune o asemenea condiție. După cum s-a văzut, aproximarea medie pătratica cere ca polinomul de aproximare sa se calculeze astfel ca valoarea abaterii

sa fie minima (fig. 3)

Fig. 3

Din punct de vedere geometric, metoda aproximării cu abatere medie pătratica minima conduce la cerința ca aria porțiunii hașurate sa fie minima.

1.4. Criteriul de aproximare în sensul lui Cebîsev

Presupunem în acest caz ca spațiul metric M este C[a,b], adică mulțimea funcțiilor continue definite pe [a,b] distanta dintre doua funcții f si g fiind data de :

Daca în locul funcție g(x) punem polinomul de aproximare Fm(x), m≤ n calculul coeficienților ck se va face din condiția ca abaterea minimax= să fie minimă. Aproximarea în acest caz se numește "în sensul lui Cebîșev". Deoarece în acest caz se impune ca un maximum de funcție sa fie minim, această cerință conduce prin alternanța celor doua cuvinte, la o alta denumire pentru aproximarea în sensul lui Cebîșev, și anume criteriul de minimax, iar abaterea respectivă se numește adeseori abatere de minimax.

Dacă valoarea minima găsita pentru abaterea de minimax este notata cu ε, atunci din punct de vedere geometric acest criteriu de aproximare impune ca graficul polinomului de aproximare sa nu depășească o banda centrata pe graficul funcției f(x), de lățime 2ε, unde ε are semnificația abaterii de minimax, (fig.4).

Fig. 4

Între cele trei criterii descrise mai sus, aproximarea prin interpolare este cea mai simpla, îinsă adesea precizia rezultatelor lasă de dorit. Criteriul de aproximare cu abatere medie pătratica minima este adesea preferat în problemele de ordin tehnic, în special daca se urmărește nu atât ordinul de mărime a erorii în fiecare punct, ci ca media erorilor sa fie cât mai mică.

Criteriul de aproximare în sensul lui Cebisev conduce în general la rezultatele cele mai precise, dar folosirea lui este deseori dificila. Se folosește în rezolvarea problemelor care impun controlul erorii în fiecare punct, de asemenea un grad sporit de precizie a rezultatelor.

În funcție de problema data se pot defini si alte criterii de aproximare, totuși cele trei descrise mai sus sunt în prezent cel mai des utilizate.

2. APROXIMAREA FUNCȚIILOR PRIN INTERPOLARE

2.1. Introducere

Interpolarea liniara a fost introdusa prima oara pentru calculul funcțiilor trigonometrice și a logaritmilor și se bazează pe aproximarea unei curbe date intre doua puncte, prin segmentul de dreapta ce unește doua puncte (fig. 5).

Fig. 5

Conform principiului aproximării prin interpolare, se impune coincidenta funcției cu polinomul de gradul întâi ales pentru aproximare, (din punct de vedere geometric, o dreapta) în nodurile de aproximare xk, x k+1.

Se scrie ecuația dreptei care trece prin punctele A[ xk, f(xk) ] si B [(Xk+1, f(xk+1)]:

Rezulta:

deci polinomul de aproximare de gradul întâi P(x) este de forma;

Reprezentând grafic funcția de mai sus obținem o dreapta numita dreapta de interpolare. Este evident ca abaterea: ∆(x) = f(x)- P(x) depinde de intervalul [ xk, xk,+1], de curbura funcției f(x) și deci de valoarea lui f"(x). Deoarece interpolarea liniara da rareori satisfacție (eroarea este mult prea mare), se folosește de obicei un polinom de grad superior (interpolare polinomiala).

2.2. Polinomul de interpolare al lui Lagrange

Fie funcția f: [x0, xn] →R data sub forma unei tabele de valori:

Punctele xk, numite și noduri de interpolare nu sunt în mod necesar egal distanțate. Se impune deci coincidenta funcției f(x) cu un polinom de aproximare de gradul n, în nodurile de interpolare, adică:

f(xk) = Pn(xk); k=0, 1, 2, ……n

Metoda lui Lagrange propune determinarea polinomului de aproximare fără a se rezolva sistemul considerat. Pentru aceasta se procedează astfel; se alege polinomul de aproximare notat Ln(x) de forma:

Ln(x) = yo φo (x) + y1 φ1(x) + …….. yn φn(x); yk = f(xk)

Unde φo (x), φ1(x),…… φn(x) , sunt polinoame de gradul n ce urmează a fi determinate impunând condițiile de interpolare.

Avem deci:

yo φo (xo) + y1 φ1(x0) + …….. yn φn(xo) = y0

yo φo (x1) + y1 φ1(x1) + …….. yn φn(x1) = y1

…………………………………………………………

yo φo (xn) + y1 φ1(xn) + …….. yn φn(xn) = yn

Rezulta:

φo (xo) = 1, φ1(x0) = 0, φk(xo) = 0, …, φn(xo) = 0

φo (x1) = 0, φ1(x1) = 1, φk(x1) = 0, …, φn(x1) = 0

…………………………………………………………………

φo (xk) = 0, φ1(xk) = 0, φk(xk) = 1, …, φn(xk) = 0

φo (xn) = 0, φ1(xn) = 0, φk(xn) = 0, …, φn(xn) = 1

Rezulta ca polinomul φk(x) de exemplu, se anulează pentru n puncte si ia valoarea 1 intr-un singur punct si anume în xk. Rezulta ca putem scrie polinomul φk(x) în modul următor;

φk(x) = ak(x – x0) (x – x1) ……(x – x k-1) (x – x k+1) ……..(x – xn)

Constanta ak se determina impunând ca în punctul xk, polinomul φk(x) sa ia valoarea 1, adică:

ak(xk – x0) (xk – x1) ……(xk – x k-1) (xk – x k+1) ……..(xk – xn) = 1,

ak = 1/(xk – x0) (xk – x1) ……(xk – x k-1) (xk – x k+1) ……..(xk – xn)

Rezultă polinomul

iar în final polinomul de interpolare al lui Lagrange devine:

Notând:

avem

Polinomul de interpolare al lui Lagrange coincide cu funcția f(x) în nodurile de interpolare, dar între două noduri poate să se abată mult de la valoarea funcției. Pentru a calcula funcția abatere se procedează astfel:

Fie funcția auxiliară:

unde k este o constantă ce urmează a fi determinată. Se observă că h(xk) = 0, pentru k = 0, 1,2…n, deci, funcția h(x) are n + 1 zerouri, oricare ar fi constanta k. Alegem această constantă impunând că h(x) să se anuleze și pentru

Rezultă

Dacă se socotește și rădăcina , rezultă că funcția h(x) are n + 2 zerouri pe (x0, xn); în plus, satisface condițiile teoremei lui Rolle pe fiecare din cele (n + 1) subintervale care se formează folosind nodurile x0, x1,…xn, . Rezultă că h’(x) are cel puțin n + 1 zerouri pe x0, xn. Raționând în mod asemănător, deducem că h”(x) are cel puțin n zerouri pe x0, xn, iar pentru derivate de ordinul n + 1 există cel puțin un punct în care aceasta se anulează

Rezultă:

Dacă egalăm cele două expresii obținute pentru valoarea lui k, rezultă:

Deoarece valoarea este oarecare în intervalul x0, xn, xk, (dar pentru = xk), evident n() = 0 se obține funcția eroare de forma

Dacă funcția f(x) are o derivată de ordinul n + 1 mărginită în valoare absolută, adică există M 0, astfel că oricare ar fi x (x0, xn), să rezulte , funcția eroare poate fi rescrisă astfel:

Deoarece în majoritatea cazurilor în care se utilizează polinomul de interpolare Lagrange, funcția f(x) este dată printr-o tabelă de valori, derivata de ordinul (n + 1) nu este cunoscută și formula anterioară este prea puțin utilă. Acesta este unul dintre dezavantajele formulei lui Lagrange, altul constând în faptul că ori de câte ori se dorește să se înlocuiască polinomul de interpolare Lagrange cu unul de grad superior, toate calculele trebuie să fie refăcute.

Aplicația 1

Să se calculeze polinomul de interpolare al lui Lagrange de gradul doi ce aproximează funcția

xk 0 1 3

yk 0 1 27

Avem:

Polinomul Lagrange de interpolare este util și pentru realizarea interpolării inverse, deoarece acest tip de interpolare folosește tabele de valori inegal distanțate.

Aplicația 2

Pentru tabela de valori dată de aplicația 1 se cere să se determine valoarea lui x corespunzătoare valorii y = 2.

Folosim polinomul de interpolare Lagrange de gradul doi, considerând valorile jR ca fiind șirul argumentelor. Avem:

Rezultă:

Polinomul Lagrange poate fi extins pentru interpolare funcțiilor de două variabile. Într-adevăr, fie (m + 1) (n + 1) puncte distincte în care este dată funcția f(xi ji), i = 0, 1, …m; j = 0, 1, … n.

Se caută un polinom L (x,y) de grad cel mult m în x și de grad cel mult n în y, astfel ca:

L(xi, yj) = f(xi, yj), i = 0, 1, … m; j = 0, 1, … n

Notând:

se demonstrează că polinomul Lagrange are în acest caz forma:

2.3. Diferențe finite

Fie funcția y = f(x) tabelată în puncte echidistante și fie h pasul tabelei.

Se numește diferența finită la dreapta pentru funcția f(x), expresia

Diferențele finite de ordin superior se definesc în mod recurent:

Se verifică imediat relațiile:

Diferențele finite de ordin superior se calculează pornind de la definiție astfel:

Se demonstrează prin inducție completă relația

Formulele de mai sus se simplifică dacă se calculează diferențele finite în punctele tabelei, deci pentru valorile x = xk.

Ținând seama că xk + h = xk+1 și notând f(xk) = yk avem:

În general:

O formulă utilă este, următoarea:

Ea este asemănătoare ca formă cu formula cunoscută din calcul integral:

Aplicația 4

Folosind formula de calcul

să se verifică că

Avem nevoie de o funcție pentru care yi = i adică o problemă similară aceleia din calculul integral în care se cere determinarea funcției y(x) astfel ca dy(x) = x.

Obținerea lui yi cunoscând pe yi poartă numele de integrare finită. Se verifică ușor că pentru y = i rezultă yi de forma yi= (i2 – i)/2. Deci

O altă formă utilă purtând numele de sumare prin părți este următoarea:

Ea corespunde formulei de integrare prin părți a calculului integral

Formula de sumare prin părți se obține observând mai întâi prin calcul direct că

Sumând de la i = 0 la i = (n – 1) și folosind formula

rezultă

Aplicația 5

Se cere să se calculeze suma seriei

Observând că se notează

Aplicând formula însumării prin părți, se obține:

Observând că ultima sumă corespunde unei progresii geometrice cu rația subunitară a, avem:

Trecând la limită, rezultă:

În aplicații calculul diferențelor finite se face mai comod prin simple operații de scădere, pornind de la tabela de valori a funcției f(x) așezată vertical. Pentru o tabelă în cinci puncte se obține:

Diferențele finite de un anumit ordin se calculează scăzând diferențele finite de ordin inferior între care se găsesc, situate la stânga.

De exemplu:

Diferențele finite: etc. se numesc diferențe finite la dreapta (progresive) . yk se calculează prin efectuarea diferențelor între ordonatele funcției y0, y1 …yn pornind de la începutul tabelei spre dreapta. Se pot defini și diferențele finite la stânga (regresive) care implică pentru calculul lor formarea diferențelor între ordonatele funcției y0, y1…, yn, pornind de la sfârșitul tabelei spre stânga. Se definesc astfel:

2.4. Polinomul lui Newton de interpolare

Dacă punctele de interpolare sunt echidistante, atunci pe baza diferențelor finite se poate construi polinomul de interpolare al lui Newton. Se vor folosi în cele ce urmează diferențele finite la dreapta. Se pune deci problema determinării unui polinom de grad mai mic sau egal cu n, care să satisfacă condițiile de interpolare.

Pn(xk) = yk, k = 0,1,2,…. n

Se ia în considerare un polinom de forma:

Impunând condițiile de interpolare polinomul se obțin coeficienții ck din ecuațiile:

Ecuațiile de mai sus pot fi puse sub o formă mai simplă dacă se ține seama că nodurile sunt echidistante.

c0 = y0

c0 + hc1 = y1

c0 + 2hc1 + 1 2h2c2 = y2

……..

c0 + nhc1 + n(n-1)h2c2 + … + 1 . 2 …nhncn = Yn

Din prima ecuație rezultă coeficientul c0 = y0. Introducând această valoare în ecuația a doua, obținem:

Procedând în mod asemănător rezultă și ceilalți coeficienți:

Obținem în final polinomul lui Newton de interpolare cu diferențe finite la dreapta.

Pornind de la un polinom de interpolare de forma:

calculele, descrise mai sus, pun în final în evidență polinomul lui Newton de interpolare cu diferențe finite la stânga:

Polinomul lui Newton cu diferențe finite la dreapta se recomandă a fi utilizat când dorim să aproximăm valorile funcției pentru x apropiat de x0 (începutul tabelei de valori), iar polinomul lui Newton cu diferențe finite la stânga va fi folosit pentru aproximarea valorilor funcției pentru x apropiat de xn (sfârșitul tabelei de valori).

Pentru a pune în evidență eroarea aproximării prin cele două polinoame de interpolare Newton să reluăm expresia erorii n(x) pentru aproximarea polinominală, pusă în evidență în paragraful precedent.

unde se găsește între x0 și xn, iar x este un punct din domeniul de definiție a lui f diferit de nodurile de interpolare. Se poate demonstra relația:

Renunțând la limită se pot face aproximările:

Din cele de mai sus rezultă expresia erorilor pentru polinoamele lui Newton de interpolare la dreapta, respectiv la stânga.

Aplicația 6

Pentru funcția dată prin următoarea tabelă de valori:

polinomul de interpolare cu diferențe finite la dreapta de gradul doi în punctul x = 1,5 dă valoarea:

De remarcat că aceeași valoare o regăsim și dacă folosim polinomul de interpolare al lui Lagrange. Eroarea care se comite prin aproximare:

Aplicația 7

Pentru funcția dată prin următoarea tabelă de valori:

polinomul de interpolare cu diferențe finite la stânga de gradul trei în punctul x = 3,5 dă valoarea:

Deoarece funcția tabelată f(x) = x3 – 2×2 + 7x – 5 este de forma unui polinom de gradul trei, valoarea calculată este valoarea exactă. De notat că tabela de valori de mai sus este identică cu tabela de valori de la aplicația 6, numai indicii valorilor xk diferă.

Se pot da și diferențe finite centrale definite astfel:

……..

Deoarece eroarea de interpolare tinde să fie mai mică atunci când se folosesc noduri de interpolare de ambele părți ale argumentului x, rezultă că diferențele finite centrale sunt indicate în special pentru interpolare în valori ale lui x situate la mijlocul tabelei de valori.

Diferențele centrale sunt de asemenea folosite cu predilecție la rezolvarea problemelor la limită.

Un tabel cu diferențe finite centrale se prezintă astfel:

Folosind în cazul tabelei de mai sus diferențele finite situate pe linia trasată continuu se obține pentru calculul polinomului de interpolare (procedând asemănător ca în cazul diferențelor finite la dreapta), formula lui Gauss la dreapta:

unde n este impar.

Dacă se utilizează diferențele finite situate pe linia trasată punctat, se obține formula lui Gauss la stânga:

unde n este impar.

Făcând media aritmetică a celor două polinoame se obține formula lui Stirling. De exemplu, folosind numai diferențele finite centrale până la ordinul 3, se obține formula:

Aplicația 8

Pentru funcția dată prin tabela de valori de mai jos să se găsească valoarea funcției în punctul x = 2,5

Folosind formula lui Gauss la dreapta, polinomul de interpolare de gradul trei în punctul x = 2,5 dă pentru funcția tabelată valoarea exactă:

2.5. Diferențe divizate. Polinoame de interpolare pe baza diferențelor divizate

Polinoamele de interpolare pot fi generalizate pentru funcții reprezentate prin tabele de valori inegal distanțate (deci fără pas constant) folosind noțiunea de diferență divizată.

Diferențele divizate de la ordinul întâi la ordinul n sunt notate și definite în felul următor:

……….

Diferențele divizate cu unul sau mai multe argumente egale pot fi definite folosind un procedeu de trecere la limită, de exemplu:

De asemenea, diferențele divizate sunt invariante față de toate permutările argumentelor xk cu condiția ca valorile yk din formula de calcul să fie permutate în același mod. Aceasta reprezintă așa – numita proprietate de simetrie a diferențelor divizate. De exemplu, deoarece:

este evident că dacă se schimbă x0 cu x1 și y0 cu y1, atunci cei doi termeni din partea dreaptă a semnului egal se schimbă între ei Deci y(x0,x1) = y(x1,x0). În cazul argumentelor xk egal distanțate, diferențele divizate se reduc la diferențe finite astfel:

Pentru calculul efectiv al diferențelor divizate se folosesc formulele de definiție și o tabela organizată după modelul tabelei de calcul a diferențelor finite.

Aplicația 9

Să se calculeze diferențele divizate pentru tabela de valori:

Avem:

Deci:

Folosind diferențele divizate, polinomul de interpolare al lui Newton se obține de forma:

Nn(x) = y0 +(x – x0) y (x0, x1) + (x – x0)(x – x1) y (x0,x1,x2) + …. + (x-x0)(x-x1) …

…(x-xn-1) y (x0, x1, …xn)

În formula de mai sus argumentele xk nu trebuie să mai fie egal distanțate. Se poate arăta cu ușurință că polinomul de interpolare al lui Newton scris mai sus, reprezintă cu adevărat un polinom de interpolare, adică:

Nn(xk) = yk, k = 0, 1, …..n. Faptul că Nn(xo) = y0 este evident. Folosind definiția diferențelor divizate și proprietatea de simetrie, putem scrie în continuare:

yk = y0 + (xk – x0)y(x0,xk)

y(x0,xk) = Y(x0, x1) + (xk – x1) y (x0, x1,xk)

y(x0,x1,xk) = y(x0,x1,x2) + (xk – x2)y(x0,x1,x2,xk)

y(x0, … xn-2,xk) = y(x0, … xn-1) + (xk – xn-1) y (x0, … xn-1, xk)

Linia a doua a fost obținută astfel:

În continuare, linia întâi poate fi folosită pentru a verifica că Nn (x1) = y1

Substituind linia a doua în linia întâi, se obține:

yk =y0 + (xk – x0)y(x0 x1) + (xk – x0)(xk – x1) y (x0, x1, x2)

Sub această formă relația de mai sus poate fi folosită pentru a verifica că Nn(x2) = y2. Se continuă în același mod pentru a verifica coincidența polinomului de interpolare al lui Newton în toate nodurile de interpolare, cu valorile yk . Rezultatele de mai sus permit folosirea pentru polinomul de interpolare a lui Newton cu diferențe divizate a formulei erorii pentru aproximarea polinomială:

Cu ajutorul diferențelor divizate, formula erorii poate fi rescrisă astfel:

Dacă se egalează coeficienții lui xk, 0 k n din cele două exprimări Lagrange respectiv Newton ale polinomului de interpolare, se obține:

Aplicația 10

Să se calculeze polinomul de interpolare de gradul doi cu diferențe divizate al lui Newton pentru tabela de valori din aplicația 9. Avem conform formulei demonstrate anterior:

N2(x) = y0 + (x – x0)y(x0,x1) + (x – x0)(x – x1)y (x0,x1,x2)

Folosind valorile diferențelor divizate calculate la aplicația 9 obținem:

3. APROXIMAREA FUNCȚIILOR CU ABATERE MEDIE PĂTRATICĂ MINIMĂ

3.1. Cazul continuu

Acest tip de aproximare folosește cu predilecție sisteme de funcții ortogonale.

Fie un sistem de funcții Pk(x), definite pe [a,b] și aparținând spațiul vectorial L2 al funcțiilor integrabile de pătrat integrabil. Sistemul este ortogonal dacă:

unde (x) 0 este funcția pondere. În particular, (x) poate fi identică cu 1. Fie funcția f(x) definită și integrabilă pe intervalul [a,b] și un polinom generalizat

Fm (x) = c0P0(x) + c1P1(x) + …. + cmPm(x) format pe baza sistemului de funcții ortogonale. Aproximarea funcției f(x) prin Fm(x) cu abatere medie pătratică conduce ca integrala

să fie minimă.

Condițiile necesare de extrem impuse funcției I de m + 1 variabile și anume C0, C1…., Cm, conduc la sistemul:

Folosind ortogonalitatea sistemului ortogonal de polinoame, sistemul de ecuații se simplifică astfel:

Fiecare ecuație implică o singură valoare pentru coeficienții Ck, deci:

Cu aceste valori ale coeficienților Ck, polinomul ce realizează aproximarea lui f(x) cu abatere medie pătratică minimă, devine:

Pentru a evalua mărimea erorii care se comite prin acest tip de aproximare, să reluăm expresia integralei I:

Al doilea termen în suma din partea dreaptă a semnului poate fi rescris

ținând seama de expresia coeficientului ck calculat. Al treilea termen în sumă se reduce, datorită ortogonalității, la valoarea

Rezultă:

Ultima relație arată că putem micșora abaterea medie pătratică

folosind polinoame generalizate de aproximare Fm(x) formate pe baza unui număr sporit de termeni. Faptul că coeficienții ck nu depind în nici un fel de indicele m are o mare importanță. Într-adevăr, aceasta permite calculul succesiv al coeficienților ck oprindu-ne la acea valoare k pentru care abaterea medie pătratică se situează sub limitele admise. Atunci când m tinde spre infinit, se demonstrează că abaterea medie pătratică tinde către zero. Dacă nu se folosesc polinoame ortogonale, atunci o schimbare a gradului polinomului generalizat de aproximare face ca toți coeficienții să fie recalculați.

Observând că I 0, ultima relație de calcul permite să se pună în evidență o inegalitate remarcabilă (inegalitatea lui Bessel).

Printre sistemele ortogonale de funcții, frecvent utilizate, sunt sistemul de polinoame Cebișev și sistemul trigonometric. Fie sistemul de polinoame Cebișev

Polinoamele Cebișev formează un sistem ortogonal de funcții în raport cu funcția pondere .

Într-adevăr:

Relațiile de mai sus se verifică cu ușurință ținând cont că:

Avem:

Se observă că I = 0 pentru m n. Dacă m = n = 0 se obține valoarea .

Dacă m = n 0, integrala devine:

Urmând teoria generală a sistemelor ortogonale rezultă că polinomul generalizat de aproximare Fm(x) format pe baza polinoamelor Cebișev:

care face integrala

să fie minimă, are coeficienții ck dați de relația:

cu excepția coeficientului c0 calculat astfel:

Fie sistemul trigonometric de funcții

Se poate arăta că acest sistem este ortogonal pe intervalul adică:

Într-adevăr, avem:

Dacă m n,

Dacă m = n,

Asemănător se pot verifica și restul relațiilor implicate de ortogonalitate. Urmărind teoria generală, rezultă că polinomul generalizat de aproximare numit și polinom trigonometric

ce aproximează funcția f(x) definită pe intervalul -N, N periodică de perioadă 2N, cu abatere medie pătratică minimă, are coeficienții dați de relațiile :

Abaterea medie pătratică devine minimă odată cu integrala I

În cazul polinomului trigonometric se obține :

Când m crește, valoarea lui I se micșorează. Se demonstrează în teoria seriei Fourier că I tinde la zero când m tinde spre infinit dacă funcția f(x) satisface condițiile lui Dirichlet.

Aplicația 1

Să se găsească polinomul trigonometric de aproximare pentru m=1 în cazul funcției y(x) = x , – x , periodică de perioadă 2. În cazul de față N = . Obținem coeficienții:

Rezultă polinomul trigonometric

3.2. Cazul discret

Fie funcția y = f(x) dată prin tabela de valori :

xk x0 x1 ……. xn

f(xk) y0 y1 ……. yn

Se cere să se aproximeze funcția printr-un polinom de gradul m n.

f(x) a0 + a1x + a2x2 + ……. + amxm = pm(x)

astfel că abaterea medie pătratică și deci suma

să fie minimă. Această metodă imaginată mai întâi de Gauss poartă numele de metoda celor mai mici pătrate. Cazul cel mai simplu îl constituie aproximarea funcției y = f(x) printr-un polinom de gradul întâi ax+b reprezentând geometric o dreaptă (dreapta de regresie). Se impune deci ca abaterea :

să fie minimă.

Condițiile necesare de extrem impuse funcției de două variabile S(a, b), conduce la:

Aceste ecuații pot fi scrise astfel :

Introducând notațiile :

se obține soluția sistemului de mai sus de forma :

Se verifică ușor că numitorul u0u2 – u12 este diferit de zero.

Să calculăm derivatele parțiale de ordinul doi pentru funcția abatere

Se observă că primele două derivate parțiale de ordinul doi sunt pozitive. În plus, deoarece

rezultă că soluția sistemului determină pentru funcția abatere S o valoare extremă și anume o valoare minimă.

Aplicația 2

Se dă funcția

xk -1 0 1 2

f(xk) 1 0 1 4

Se cere dreapta de regresie ce aproximează funcția de mai sus.

Sistemul de ecuații ce determină coeficienții a și b ai dreptei de regresie devine în acest caz : 4b + 2a = 6, 2b + 6a = 8. Se obține soluția sistemului a = 1, b = 1, deci dreapta de regresie cerută este y = x+1.

Aplicația 4

Să se aproximeze funcția f dată prin tabela de valori calculată în 5 puncte :

xk-2 xk-1 xk xk+1 xk+2

yk-2 yk-1 yk yk+1 yk+2

conform criteriului celor mai mici pătrate, apoi să se folosească polinomul de

aproximare găsit pentru a determina derivata funcției f în punctul xk. Să facem schimbarea de variabilă (x – xk)/h = t, adică x = xk + th. Valorile x din tabela de valori se obțin pentru t = -2, -1, 0, 1 și 2. Aproximând pe f prin polinomul a0 + a1t + a2t2 se ajunge conform criteriului celor mai mici pătrate la sistemul :

Rezolvând sistemul de ecuații de mai sus se obține :

a0 = 1/35(-3yk-2 + 12yk-1 + 17 yk + 12yk+1 – 3yk+2)

a1 = 1/10(-2yk-2 – yk-1 + yk+1 + 2yk+2)

a2 = 1/14(2yk-2 – yk-1 – 2yk – yk+1 + 2yk+2)

Deoarece

f'(xk) P'(xk) = a1/h

se obține în final

Sistemul de ecuații care conduce la determinarea coeficienților polinomului de aproximare se rezolvă cu ușurință pentru valori mici ale lui m, așa cum s-a văzut, dar pentru valori mari ale lui m devine adeseori nestabil.

Această dificultate poate fi evitată lucrând cu sisteme de funcții ortogonale.

Se poate arăta astfel că dacă se folosește o tabelă de valori

x 0 1 …… 2N

y(x) y(0) y(1) ….. y(2N)

calculată într-un număr impar de puncte 2N + 1, egal distanțate, atunci polinomul trigonometric de aproximare (prin metoda celor mai mici pătrate) pentru funcția y(x) dată tabelar și periodică de perioada 2N este de forma

coeficienții ak și bk fiind calculați cu formulele :

, k = 0, 1, 2, ……….. N

, k = 0, 1, 2, ………….. N

Într-adevăr, să luăm un polinom trigonometric de aproximare oarecare, dar având forma :

N

și să arătăm că abaterea

devine minimă atunci când coeficienții Ak = ak și Bk = bk

Termenul general al abaterii este de forma :

Ridicând la pătrat și folosind proprietatea de ortogonalitate a sistemului finit de funcții trigonometrice

se obține :

Evident pentru a micșora abaterea S, trebuie să luăm A0 = a0, Ak = ak, Bk = bk

Se poate arăta că abaterea S descrește atunci când M crește atingând valoarea 0 pentru M = N, cu alte cuvinte în acest ultim caz polinomul trigonometric de aproximare coincide cu un polinom generalizat de interpolare. De multe ori se preferă să se folosească pentru funcția y(x) periodică de perioadă 2N, o tabelă de valori de forma

x 0 1 2N – 1

y(x) y(0) y(1) y(2N – 1)

dată într-un număr par de puncte 2N. Atunci polinomul trigonometric este de forma :

unde coeficienții ak, bk se calculează cu formulele :

Aplicația 5

Să se determine un polinom trigonometric de aproximare pentru funcția y(x) dată prin tabela de valori

x 0 1 2 3

y(x) 0 1 4 9

Avem 2N – 1 = 3, deci N = 2

Polinomul trigonometric va fi de forma :

Coeficienții ak, bk calculați din formulele de mai sus sunt :

a0 = 7 ; a1 = -2 ; b1 = -4 ; a2 = -3 ;

Rezultă :

Frecvent intervin în tehnică funcții periodice de perioadă 2 definite printr-o tabelă de valori dată în n puncte ce împart intervalul [0, 2 în n – 1 părți egale.

k 0 = 1, 2 …….k = (k – 1) 2/(n – 1), …….n = 2

(k) 1, 2 ……..k ……..…. n

De fapt cazul în care o funcție g(x) este definită pe intervalul , P poate fi redus la cazul precedent, printr-o schimbare de variabilă de forma :

Într-adevăr, această substituție transformă funcția într-o funcție () definită pe intervalul 0, 2. Dacă tabela este dată într-un număr par de puncte, adică n = 2N atunci polinomul trigonometric de aproximare pus în evidență în acest caz, apare sub o formă simplificată :

Coeficienții ak și bk se calculează cu formulele :

k = 0, 1, 2 ……. N

k = 1, 2 ……. N – 1

Dacă funcția f() se dă printr-o tabelă de numai 12 valori, adică n = 12, apar mari simplificări în calculul coeficienților Fourier cu relațiile scrise anterior. Într-adevăr, în calcule nu mai figurează decât 24 de funcții trigonometrice care cu excepția semnului, nu au decât trei valori distincte diferite de zero : 0,5 /2 = 0,866, 1 .

În acest caz calculul coeficienților Fourier este indicat să se facă prin metoda Runge folosind tabele tip. Calculul se va desfășura în felul următor : se scriu cele 12 valori fi ale funcției f astfel :

Sumele Sk se calculează în modul următor : S0 = f12, S6 = f6, în rest Sk reprezintă suma valorilor funcție f aflate pe coloana k. Diferențele dK se formează luând de asemenea în considerare valorile funcției aflate pe coloana K în tabelă.

De exemplu d3 = f3 – f9. Se calculează în continuare în mod analog sumele k și k respectiv diferențele k și k folosind tabelele.

Calculul coeficienților Fourier se face în continuare utilizând următoarele două tabele:

Coeficienții sinusurilor

Coeficienții cosinusurilor

În cele două tabele cifrele 0,500, 0,866 și 1,000 indică faptul că trebuie multiplicate prin acestea valorile corespunzătoare aflate pe linia ce corespunde cifrei respective în coloanele din partea dreaptă. se calculează în continuare sumele I și II ca în tabelele de sume. În final rezultă coeficienții ak și bk luând în considerare valorile respective I + II, de. asemenea I – II. Se pot scrie de exemplu, folosind tabelele de mai sus :

6b1 = 0,5 . 1 + 1 . 3 + 0,866 2

6b3 = 1 . 1 – 1 . 3

6a2 = – 0,5 2 – 1 . 3 + 0,5 1 + 1 . 3

Pornind de la sistemul de polinoame Cebîșev {Tn (x)} se poate forma polinomul de aproximare pentru funcția y (x) dată prin tabela de valori :

unde xk = [(2k + 1) / (2N)] k = 0,1 ………… N -1

de forma : P(x) = C0 T0 (x) + c1 T1 (x) + ……….+ cm Tm (x).

Se observă că valorile xk sunt alese ca fiind zerourile polinomului Cebișev TN (x). Coeficienții ck se calculează folosind proprietatea de ortogonalitate pentru sistemele finite de polinoame Cebișev.

0, i j

N/2, i = j 0

N, i = j = 0

Coeficienții se obțin de forma :

c0 = ck =

4. APROXIMAREA FUNCȚIILOR ÎN SENSUL LUI CEBÎȘEV

4.1. Cazul continuu

Polinomul de gradul n poate fi considerat un polinom generalizat format pe baza monoamelor 1, x, x2, …… xn. Examinând comportarea acestor monoame pe intervalul -1, 1, se observă că fiecare din ele își atinge valoarea maximă 1, la x = ±1 și valoarea minimă 0, la x = 0. Cu alte cuvinte valorile extreme sunt atinse la capetele și la mijlocul intervalului. Dacă funcția f(x) este aproximativă prin polinomul : f(x) a a1x a2x2 ….. anxn, modificarea coeficienților sau schimbarea gradului polinomului, operații care se execută în procesul de aproximare, vor produce erori mici pentru valori ale lui ''x'' situate în apropierea lui zero și erori mari pentru valori ale lui ''x'' situat în apropierea lui ±1. Deci, cu toate avantajele care provin din expresia lor simplă, setul de monoame are dezavantajul unor extreme neuniform distribuite. Mult mai potrivite par din acest punct de vedere funcțiile trigonometrice : cos , cos 2, ………. cos n, așa cum se observă din fig. 6.

Faptul că aceste funcții sunt definite pe intervalul 0, nu constituie o piedică, știut fiind că putem transforma un interval oarecare a, b în intervalul -1, 1 prin schimbarea de variabilă

Fig. 6.

Dificultatea utilizării funcțiilor cos n în aproximarea funcțiilor constă în faptul că sunt dificil de evaluat numeric. Din această cauză se preferă transformarea funcțiilor cos n definite pe intervalul -, într-un set de polinoame definite pe intervalul -1, 1, folosind transformarea : = arc cos x

Obținem : T1(x) = cos = cos(arc cos x) = x

T2(x) = cos 2 = cos(2arc cos x) = 2 cos2 (arc cos x) – 1 = 2×2 – 1

Tn(x) = cos n

Pentru calculul polinoamelor Tn(x), numite polinoame Cebîșev, se folosește identitatea trigonometrică

cos n = 2 cos x cos (n – 1) – cos (n – 2)

care conduce la relația de recurență :

Tn(x) = 2x . Tn-1(x) – Tn-2(x)

Primele opt polinoame Cebîșev sunt următoarele :

T0 = 1 T4 = 8×4 – 8×2 1

T1 = x T5 = 16×5 – 20×3 5x

T2 = 2×2 – 1 T6 = 32×6 – 48×4 18×2 – 1

T3 = 4×3 – 3x T7 = 64×7 – 112×5 56×3 – 7x

Din relațiile de mai sus rezultă cu ușurință puterile lui x în funcție de polinoamele Cebîșev

1 = T0

x = T1

x2 =

x3 =

Folosind definiția, rezultă că Tn(x) este zero pentru = (2k + 1)/2n, deci rădăcinile polinomului sunt :

Cele n valori de mai sus sunt cuprinse între -1 și 1. Ținând seama că Tn(x) are numai n zerouri, rezultă că toate sunt cuprinse în intervalul -1, 1. Deoarece polinomul Tn(x) este egal cu un cosinus, rezultă că în intervalul -1, 1 nu poate să depășească valoarea 1, atingând aceste valori limită în (n + 1) puncte.

Tn(x) = (-1)k, în punctele k = cos k/n ; k = 0, ……. n

Ansamblul punctelor k formează suportul lui Cebîșev. O proprietate utilă a polinoamelor Cebîșev o constituie faptul că dintre toate polinoamele Pn(x) de gradul n care au coeficientul lui xn egal cu 1, cel pentru care abaterea sa de la zero este minimă în intervalul -1, 1 este polinomul

21-nTn(x)

Într-adevăr, să observăm că deoarece valoarea maximă absolută pentru Tn(x) este 1, atunci pentru 21-nTn(x) această valoare este 21-n ; în plus, polinomul 21-nTn(x) are coeficientul termenului xn egal cu 1. Să presupunem acum că ar exista un polinom pn(x) având coeficientul lui xn egal cu 1, astfel ca pn(x) 21-n, pentru orice x -1, 1. Atunci datorită faptului că polinomul 21-nTn(x) ia valori extreme alternante ca semn în cele n puncte, k = cos k/n, avem :

pn(0) 21-nTn(0), pn(1) 21-nTn(1), etc.

Rezultă că polinomul diferență pn(x) – 21-nTn(x) are gradul cel mult n-1, deoarece coeficienții lui xn în cele două polinoame pn(x) și 21-nTn(x) sunt egali, în plus polinomul diferență schimbă semnul de n ori în fiecare din cele n intervale (k-1, k), k = 0,1,2 … n-1. Se ajunge deci la o contradicție. Rezultă că nu există un polinom p(x) având abaterea sa de la zero mai mică decât 21-nTn(x). Cele de mai sus pot fi utilizate pentru micșorarea utilizării erorii ce apare în aproximarea unei funcții printr-un polinom de interpolare de gradul n. În acest caz, eroarea este :

Singura cale efectivă de a micșora eroarea este de a minimiza valoarea maximă a polinomului de gradul (n+1) care apare în partea dreaptă a relației anterioare. Acest lucru implică alegerea nodurilor de interpolare x0, x1, ……. xn ca fiind zerourile polinomului Cebîșev Tn+1(x), adică :

Dacă se lucrează pe un interval arbitrar ayb, deci se trece de la intervalul -1, 1 la intervalul a, b, nodurile de interpolare yk se obțin din valorile xk prin relația de transformare

Teoria polinoamelor Cebîșev poate fi generalizată. Se demonstrează astfel că dacă funcția f(x) este continuă pe intervalul a, b, există un polinom P(x) unic, de gradul n numit polinom de minimax, astfel că abaterea :

minimax = max f(x) – Pn(x)

xa, b

este minimă. În plus, există cel puțin (n+2) puncte în care minimax atinge valoarea sa maximă, alternând semnul acesteia. Această observație deschide și calea pentru calculul polinomului de minimax de gradul n.

Pn(x) = c0 + c1x + c2x2 + ….. + cnxn

Într-adevăr, fie x0, x1 ….. xn+1 cele (n+2) puncte în care funcția abatere atinge valorile sale limită ±M. Se poate scrie sistemul de n+2 ecuații :

f(x0) – Pn(x0) = M

f(x1) – Pn(x1) = -M

f(xn+1) – Pn(xn+1) = (-1)nM, unde = +1 sau -1

În plus, deoarece valorile x0, x1, x2, …… xn corespund punctelor în care funcția abatere ia valoare extremă, derivata sa în aceste puncte trebuie să fie egală cu zero.

f '(x0) – P'n(x0) = 0

f '(x1) – P'n(x1) = 0

f '(xn+1) – P'n(xn+1) = 0

Punând împreună relațiile de mai sus, obținem un sistem de 2n+4 ecuații ce poate fi scris prescurtat astfel :

f(xk) – Pn(xk) = (-1)kM

f '(xk) – P'n(xk) = 0, k = 0,1, …..n+1

= +1 sau -1

Sistemul celor 2n+4 ecuații conțin 2n+4 necunoscute și anume :

c0, c1, ……….., cn, x0, x1, ……….., xn+1, M

Aplicația 1

Să se determine polinomul de minimax de gradul 1 care aproximează funcția f(x) cu f''(x) 0, pe intervalul a, b. Să presupunem polinomul de minimax căutat, de forma P(x) = mx + n. Conform teoriei, trebuie să găsim trei puncte x1 x2 x3 în a, b pentru care funcția abatere (x) = f(x) – P(x) atinge valorile sale extreme alternând ca semn. Deoarece '(x2) = f'(x2) – P'(x2) = f'(x2) – m = 7, rezultă f'(x2) = m. Dar f''(x) 0, deci f' este strict crescătoare și poate să fie egală cu m doar o singură dată. Rezultă că în interiorul intervalului se găsește doar punctul x2 ; x1 și x3 sunt situate la capetele intervalului, x1 = a, x3 = b.

Rezultă sistemul de ecuații :

f(a) – P(a) = -f(x2) – P(x2) = f(b) – P(b) = M,

unde M reprezintă abaterea maximă.

Rezolvând acest sistem obținem :

unde punctul x2 este determinat din relația

De exemplu pentru funcția f(x) = x2, x 0, 2 se obține f’(x2) = 2×2 = 2, deci x2 = 1

În continuare, folosind formulele de mai sus rezultă m = 2, n = -1/2.

Fig. 7.

În concluzie, polinomul de minimax de gradul întâi este P(x) = 2x – 1/2, având reprezentarea grafică dată în figura 7. Se observă că graficul polinomului de minimax nu depășește o bandă trasată centrat pe graficul funcției de lățime egală cu 2M = 1. Valoarea maximă a abaterii M = -1/2 se obține în punctul x = 0.

Se poate compara pe fig. 7. calitatea aproximării prin dreapta de minimax în comparație cu dreapta de interpolare y = 2x trasată în punctele de coordonate (0,0) și (2,4) sau cu dreapta de regresie care în cazul acesta este y = 2x – 1/3.

Dacă funcția f(x) este continuă, rezultă și minimax de asemenea funcție continuă și deoarece ia valori egale și opuse ca semn în n+2 puncte, se va anula de asemenea în n+2 puncte. Rezultă că pe intervalul (a,b) există n+1 puncte în care funcția f(x) coincide cu polinomul de minimax, deci acest polinom este în realitate un polinom de interpolare de gradul n, nodurile de interpolare fiind situate între punctele de extrem. În acest caz ne putem folosi de formula erorii pentru interpolare. Dacă derivata f(n+1)(x) există și este constantă, polinomul de minimax de gradul n poate fi construit ca un polinom de interpolare folosind, o bază de noduri de interpolare dată de zerourile polinomului Cebîșev de gradul n. Calculul exact al polinomului de minimax de grad superior lui unu este dificil.

Se poate însă arăta că dacă funcția f : -1, 1R poate fi dezvoltată într-o serie formată pe baza polinoamelor Cebîșev

,

atunci o sumă parțială.

reprezintă o bună aproximare pentru polinomul de minimax de gradul n. Dacă funcția f(x) care trebuie să fie aproximată este ea însăși de forma unui polinom, atunci în suma parțială pn(x) formată pe baza polinoamelor Cebîșev trebuie să renunțăm la ultimul termen.

Aplicația 2

Să se aproximeze funcția f(x) = cos x pe intervalul -1, 1 printr-un polinom de minimax de gradul 1.

Deoarece putem dezvolta funcția cos y în serie de puteri

avem folosind tabela care echivalează puterile lui x în raport cu polinoamele lui Cebîșev.

Neglijând în continuare pe T4 care contribuie în mică măsură la formarea sumei, avem :

Înlocuind polinoamele Cebîșev T0 și T2 prin polinoamele în x corespunzătoare, obținem expresia aproximativă a polinomului de minimax de gradul doi care aproximează funcția cos x :

În punctul x = 0 de exemplu, se obține o valoare aproximativă folosind polinomul astfel calculat egală cu 0,99479, față de valoarea exactă 1, deci o eroare egală cu 0,00521.

Aplicația 3

Să se aproximeze funcția F(y) = y2 cu un polinom de minimax pe intervalul 0, 1. Se folosește mai întâi transformarea

x = 2y – 1, y 0, 1, x -1, 1

pentru a transforma funcția F(y) definită pe intervalul 0, 1 într-o nouă funcție f(x) = (1/4)(x2 + 2x + 1) definită pe intervalul -1, 1.

Folosind tabela care face legătura între monoamele xn și expresia corespunzătoare scrisă în funcție de polinoamele Cebîșev, obținem :

Neglijând, conform teoriei, ultimul termen, rezultă :

f(x)

Revenind la intervalul [0,1] se obține o aproximare a polinomului de minimax de forma: F(y) y – 1/8

4.2. Cazul discret

Fie funcția f(x) dată într-o tabelă de valori calculate în trei puncte diferite

Se poate arăta că există o dreaptă y = mx + n care aproximează funcția dată astfel că în cele trei puncte funcția abatere ia valori egale ce alternează ca semn.

(x1) = y1 – f1 = h

(x2) = y2 – f2 = – h

(x3) = y3 – f3 = h

Într-adevăr, deoarece în punctele xk, ordonatele dreptei sunt yk = mxk + n, se verifică egalitatea:

(x3 – x2) y1 – (x3 – x1)y2 + (x2 – x1)y3 = 0

În relația de mai sus să considerăm :

y1 = f1 + h, y2 = f3 – h, y3 = f3 + h

Cu alte cuvinte, se alege dreapta ce trece prin punctele :

(x1 . f1 + h), (x2 , f2 – h), (x3, f3 + h),

deci are proprietatea menționată mai sus referitoare la funcția abatere. Rezultă :

(x3 – x2) (f1 + h) – (x3 – x1) (f2 – h) + (x2 – x1) (f3 + h) = 0

Se obține :

h =

De fapt pentru construcția dreptei sunt necesare doar două puncte, dar se verifică ușor că și cel de-al treilea punct este situat pe aceeași dreaptă. Dreapta astfel calculată poartă numele de dreapta lui Cebișev. Se verifică că dreapta lui Cebișev este unică, în plus este identică cu dreapta minimax. Aceasta înseamnă că dreapta aproximează funcția dată în sensul lui Cebișev, deci nu există o altă dreaptă care să realizeze o valoare maximă a abaterii mai mică, alternând ca semn în cele trei puncte. În cazul în care funcția este dată printr-o tabelă în mai mult de trei puncte, pentru calculul dreptei de minimax se recomandă o metodă iterativă. Aceeași metodă se folosește și în cazul căutării polinomului de minimax de grad superior.

Se ilustrează metoda de calcul iterativ pe următorul exemplu :

Aplicația 4

Se dă funcția f(x) = IxI printr-o tabelă de valori calculată în 5 puncte :

Se cere să se determine polinomul de minimax de gradul doi care aproximează această funcție folosind metoda iterativă.

PASUL 1

Metoda începe prin alegerea unui set inițial de patru puncte dintre cele date, de exemplu -2, -1, 0,1.

PASUL 2

Se calculează în punctele selectate polinomul de gradul doi p2(x) = a + bx + cx2, conform teoriei lui Cebișev deci astfel ca abaterile (xk) = p2(xk) – fk să aibă o valoare maximă egală (fie aceasta k) alternând ca semn.

Avem

a – 2b + 4c – 2 = h

a – b + c – 1 = – h

a – 0 = h

a + b + c -1 = – h

Soluția sistemului de ecuații va fi a = 1/4, b = 0, c = 1/2, h = 1/4.

Deci polinomul p2 (x) = 1/4 + x2/2.

PASUL 3

Se calculează abaterile în toate cele cinci puncte, fie acestea mk . Se notează prin M valoarea maximă absolută a acestora. Dacă IhI = M procedura se termină, polinomul de aproximare fiind acela calculat la pasul 2. Este cazul exemplului de față, în care cele cinci abateri calculate sunt: 1/4, – 1/4, 1/4, -1/4, 1/4, deci în care valoare maximă absolută este egală cu valoarea h calculată și anume 1/4. Dacă acest fapt nu ar fi avut loc, algoritmul ar fi trebuit continuat în felul următor:

PASUL 4

Se alege un nou set de patru puncte, deci se modifică setul inițial în felul următor : se adaugă un nou punct la setul inițial și anume acela în care eroarea calculată este în valoare absolută, maximă. Apoi se renunță la un punct din setul inițial și anume astfel încât noul set de patru puncte să prezinte erori având semne alternante.

Sunt reluate pasurile 2 și 3. Se poate demonstra convergența procedurii descrise mai sus către polinomul de minimax.

În cazul exemplului numeric prezent, parabola de minimax calculată aproximează funcția modul dată tabelar conform fig. 8.

Fig. 8.

Prezentăm un program pentru determinarea polinomului de interpolare în forma Newton.

include stdio.h

include conio.h

/ calculează diferența divizată f[x0, x1,…, xk]/

double diferența_divizată (int k, double x, double y);

void main (void)

double x;

double y [100];

int i, j,n;

printf(“\nn=”);scanf(“%d”, &n);

for (i = 0; i <=n; i ++)

printf(“\n x[%i]=”,i);scanf(“%lf”, &x[i]);

for (i=0; i<=n; i++)

printf(“\n y[%i]=”,i);scanf(“%lf”, &y[i]);

printf(“Poliomul Newton asociat punctelor date: “);

for (i = 0, i <=n; i ++)

printf(“\n%6.3lf”, diferența_divizată (i, x, y));

for (j = 0, j < i, j++)

printf(“(x+%6.3lf)”, x[j]);

if (i<n)

printf(“+”);

getch ();

double diferența_divizată (int k, double *x, double *y)

int i, j;

double s,p;

s=0;

for (i = 0, i <=k; i ++)

p=1;

for (j =0, j<=k; j ++)

if (j! = i)

p=p*(*(x+i) – (*(x+j)));

s+ = y[i]/p;

return s;

Fie integrabilă pe și

Dacă se cunoaște o primitivă F a lui f atunci, după cum se știe, formula Leibniz – Newton J = F(b) – F(a) ne dă valoarea exactă a integralei. Dacă însă, fie f nu admite primitive, fie f nu admite o primitivă simplă, fie f nu admite primitive elementare, atunci se aproximează f cu o altă funcție g așa încât integrala acesteia să se poată calcula exact și luăm

În cele ce urmează vom ilustra această metodă luând ca funcție g de aproximare polinomul Lagrange asociat funcției date f și considerând x0 = a, xn = b.

Avem

și notând obținem formula aproximativă

(1)

numită formula lui R.Côtes. Remarcăm faptul că numerele ci nu depind de funcția f și pot fi tabelate.

Folosind teorema 2, § 1, se poate da și o estimare a erorii absolute în formula aproximativă Côtes.

rezultă deci

(2)

unde

Cazuri particulare ale formulei Côtes

1) Formula dreptunghiurilor. Vom aproxima f C1 printr-un polinom Lagrange de grad 0, deci pentru n = 0, formula Côtes )1= devine

adică

(3)

care se numește formula dreptunghiului, întrucât aria subgraficului funcției f este aproximată cu aria unui dreptunghi, figura 9.

Din (2) obținem

și deci

(4)

unde

Fig. 9.

Să împărțim intervalul în n părți egale și să aplicăm formula dreptunghiului (3) fiecărui interval Obținem

adică

(5)

care reprezintă formula dreptunghiurilor.

Folosind (4) deducem eroarea absolută

Fig. 10.

și cum este pasul diviziunii, rezultă că eroarea în formula (5) este de ordinul h2.

, unde (6)

Dacă se cere calculul integralei cu o eroare mai mică decât 0, atunci punem condiția

de unde rezultă

(6’)

2) Formula trapezelor. Vom aprecia funcția f C2 printr-un polinom Lagrange de grad unu, deci pentru n = 1, formula Côtes (1) devine

în care y0 = f(a),

deci obținem

(7)

care se numește formula trapezului, întrucât aria subgraficului funcției f se aproximează cu aria unui trapez, figura 6.

Fig. 11.

Din (2) obținem

adică

, unde (8)

Împărțind intervalul în n părți egale aplicând formula trapezului (7) fiecărui subinterval , i = 0, n – 1 obținem

Fig. 12.

adică

(9)

care reprezintă formula trapezelor

Din (8) deducem eroarea absolută în formula generală a trapezelor

adică

, unde (10)

Dacă se cere calculul integralei cu o eroare mai mică decât 0, atunci punem condiția

de unde rezultă

(10’)

3. Formula Simpson (formula parabolelor). Vom aproxima funcția f C3cu polinomul Lagrange de grad 2, luând nodurile deci pentru n = 2, formula Côtes (1) devine

(11)

unde

și deci (11) devine:

(12)

numită formula celor trei nivele sau formula parabolei, întrucât, aria subgraficului funcției f este aproximată prin aria subgraficului parabolei care trece prin punctele ca în figura 13.

Fig. 13.

Din (2) obținem eroarea absolută

,

adică

, unde (13)

Să împărțim acum intervalul în n = 2 m părți egale și să aplicăm formula celor trei nivele (12) fiecărui interval obținem

adică (14)

care reprezintă formula generală Simpson sau formula parabolelor

Din 13 obținem eroarea absolută a formulei Simpson (14),

deci

unde (15)

Dacă se cere calculul integralei cu o eroare mai mică decât 0, atunci punem condiția

de unde rezultă

(15)

Observație

Dacă atunci se poate obține următoarea evaluare a erorii în formula (14)

unde

Exemple

1. Luând n = 4, să se calculeze cu aproximație integrala folosind formulele: dreptunghiurilor, trapezelor, parabolelor (Simpson).

Formula dreptunghiurilor (5) pentru n = 4 devine

Avem , și deci, conform (6), eroarea este mai mică decât

Formula trapezelor (9) în acest caz devine

Avem și deci, conform (10), eroarea este mai mică decât

c) Formula Simpson (14) în cazul de față devine

Avem:

și deci, conform (16), eroarea este mai mică decât .

Se observă că din cele tei formule utilizate, a dreptunghiurilor, a trapezelor și a lui Simpson, la calculul integralei , luând n = 4, cel mai bun rezultat a fost obținut cu formula lui Simpson, eroarea în acest caz fiind mai mică decât 10-3.

Întrucât rezultă ln 2 0,69324 cu o eroare mai mică decât 10-3.

20. Să se calculeze cu o eroare mai mică de 10-2 integrala folosind metoda trapezelor.

Avem:

deci

Din tabloul de variație a funcției f” rezultă că f” are un minim egal cu – 2 în 0 și un maxim egal cu în 1.

Deci

Determinăm n 1 din condiția care ne dă n 5.

Atunci avem, conform (9)

Deci cu eroare mai mică decât 10-2 și cum , obținem = 3,13492 cu o eroare mai mică decât .

Având în vedere complexitatea calculelor numerice care apare în metodele prezentate, dăm câte un program în C pentru fiecare din cele trei formule (dreptunghiurilor, trapezelor și respectiv Simpson) pe care le ilustrăm pe o integrală concretă.

Exemplu

Să se calculeze cu o eroare mai mică decât = 10-2, folosind fiecare din cele trei metode, integrala:

Avem:

Derivatele fiind pozitive rezultă că f’, f’’, f’’’ sunt crescătoare, deci avem:

a) În formula dreptunghiurilor numărul natural n se determină din (6’) care, în acest caz devine

și putem lua n = 272

b) În formula trapezelor numărul n se determină cu (10’), adică

și putem lua n = 12

c) În formula Simpson, numărul n se determină din (15’), adică

și putem lua n = 4.

Prezentăm acum cele trei programe și valoarea integralei astfel obținută.

Program pentru metoda dreptunghiurilor:

include stdio.h

include math.h

double f (double x)

{

return exp (x* x);

}

void metoda _ dreptunghiurilor (double a, double b, double M, double epsilon)

int n, i;

double s, x, pas;

n = floor (M* (b – a) * (b – a) / (2 * epsilon)) + 1;

s = 0;

x = a;

pas = (b – a) / n;

for (i = 0; i n; i ++)

s + = f(x);

x + = pas;

print (“\ nn = % d\n Valoarea integralei este: % 6.3 lf”, n,s* pas);

void main ()

double a,b, M, epsilon;

/* citesc a,b,M, epsilon */

printf (“\na = “);

scanf (“%lf”, &a);

printf (“\nb = “);

scanf (“%lf”, &b);

printf (“\nM = “);

scanf (“%lf”, &M);

printf (“\nepsilon = “);

scanf (“%lf”, & epsilon);

metoda_dreptunghiurilor (a, b, M, epsilon);

Rezultatul dat de program este I 1,459

Program pentru metoda trapezelor:

include stdio.h

include math.h

double f (double x)

{

return exp (x* x);

}

void metoda _ trapezelor (double a, double b, double M, double epsilon)

int n, i;

double s, x, pas;

s = b – a;

n = floor (sqrt (M* s* s* s/ ( 12* epsilon))) + 1;

s = 0;

pas = (b – a) / n;

x = a + pas;

for (i = 1; i n; i ++)

s + = f(x);

x + = pas;

print (“\ nn = % d\n Valoarea integralei este: % 6.3 lf”,n, (f(a) + f(b) + 2*s)* pas/2);

void main ()

{

double a,b,M epsilon;

/*citesc a,b,M epsilon */

printf (“\na = “);

scanf (“%lf”, a);

printf (“\nb = “);

scanf (“%lf”, &b);

printf (“\nM = “);

scanf (“%lf”, &M);

printf (“\nepsilon = “);

scanf (“%lf”, & epsilon);

metoda _trapezelor (a, b, M, epsilon);

Rezultatul dat de program este I 1,466.

c) Programul pentru metoda Simpson:

include stdio.h

include math.h

double f (double x)

{

return exp (x* x);

}

void metoda _ Simpson (double a, double b, double M, double epsilon)

int n, i;

double s1, s2, x, pas, pas2, t;

t = b – a;

t = t * t;

n = floor (pow(M* t* t* / ( 192* epsilon), (double) 1/3)) + 1;

s1 = 0;

s2 = 0;

pas = (b – a) / n;

pas2 = pas + pas;

x = a + pas;

for (i = 1; i n; i + = 2)

s1 + = f(x);

x + = pas2;

x = a + pas 2;

for (i = 2; i n; i + = 2)

s2 + = f(x);

x + = pas2;

print (“\ nn = % d\n Valoarea integralei este: % 6.3 lf”,n, (f(a) + f(b) + 4*s1 + 2* s2)* pas/3);

void main ()

{

double a,b,M epsilon;

/*citesc a,b,M epsilon */

printf (“\na = “);

scanf (“%lf”, a);

printf (“\nb = “);

scanf (“%lf”, &b);

printf (“\nM = “);

scanf (“%lf”, &M);

printf (“\nepsilon = “);

scanf (“%lf”, & epsilon);

metoda _Simpson (a, b, M, epsilon);

Rezultatul dat de program este I 1,464.

BIBLIOGRAFIE

Larionescu, D., Metode numerice și programarea în limbajul FORTRAN, Institutul Politehnic București, 1985.

Nicolae, P., Otlăcan, P. Elemente de analiză numerică și programarea calculatoarelor electronice. București, Academia Militară, 1971

Alexandru Petcu, Horia Cristian Petcu, Analiza Numerică, Editura Premier, Ploiești, 2001

Alexandru Petcu, Analiza Numerică, Editura Universității Petrol-Gaze Ploiești, 1997

Stănășilă, O. Analiză matematică. București, Ed. Didactică și Pedagogică, 1981.

Șabac, I. Gh. ș.a. Matematici speciale. București, Ed. Didactică și Pedagogică, 1983.