APROXIMAREA POLINOMIALA PRIN INTERMEDIUL FUNCTIILOR SPLINE CUBICE NATURALE

APROXIMAREA POLINOMIALA PRIN INTERMEDIUL FUNCTIILOR SPLINE CUBICE NATURALE

In foarte multe aplicatii apare necesitatea evaluarii aproximative unei functii . Se considera o diviziune a intervalului [a,b]:

pentru care se cunosc valorile functiei :

Se pune problema determinarii unui polinom, P, de grad maxim (n - 1), care sa treaca prin cele n puncte (a se revedea notele de curs):

In prezenta lucrare se abordeaza interpolarea polinomiala cu functii spline cubice naturale. Aceasta va fi folosita si pentru evaluarea primelor doua derivate ale unei functii cunoscuta printr-un sir de valori, precum si pentru evaluarea integralei definite pe un interval finit.

1. Interpolarea si aproximarea numerica cu ajutorul functiilor spline cubice naturale

Functiile spline cubice naturale sunt functii de tip polinomial, definte pe ramuri prin intermediul unor polinoame de gradul al III-lea pentru fiecare interval :

(1)

unde seturile de coeficienti se determina din conditiile:

(C1)          conditii de interpolare

(C2)          continuitatea functiei :

(C3)          continuitatea primei derivate, :

(C4)          continuitatea derivatei a doua, :

(C5.1) ;

(C5.2) .

Pentru , , valoarea calculata aproximeaza valoarea cautata a functiei f: . Prin derivarea relatiei (1) se pot obtine aproximatii pentru derivatele lui f evaluate in , iar prin integrarea lui (1) pe un interval [qa, qb], a ≤ qa < qb ≤ b (cu notatiile anterioare din prezenta lucrare), se obtine valoarea aproximativa a integralei definite a functiei f pe intervalul dorit. Rezulta urmatoarele formule de calcul, pentru fiecare interval :

pentru aproximarea derivatei I-a :

(2)

pentru aproximarea derivatei a II-a:

(3)

pentru aproximarea integralei definite pe intervalul :

(4)

Se prezinta in continuare un grup de functii MATLAB care permit implementarea conceptelor expuse.

2. Functii MATLAB pentru implementarea aproximarii numerice prin intermediul functiilor spline cubice naturale

Functia MATLAB splinen

Scop Interpolarea folosind functii spline cubice naturale

Apel ye = splinen(xi,yi,xe)

sau

pp = splinen(xi,yi)

Parametri formali

Ø INTRARE : (xi,yi) - vectori continand punctele (argument, valoare functie in argument) pentru care se doreste interpolarea;

- valorile din vectorul xi trebuie sa fie distincte, altfel rutina este abandonata;

- trebuie sa existe cel putin patru valori in fiecare din vectorii xi si yi, altfel rutina fiind abandonata .

xe - vectorul de argumente pentru care se doreste aproximarea functiei (numai pentru prima forma de apel).

Ø IESIRE : ye - valorile aproximante calculate, corespunzatoare valorilor argumentului, furnizate in vectorul de intrare xe (numai pentru prima forma de apel).

pp - matrice continand argumentele (nodurile) xi, precum si coeficientii de interpolare calculati (numai pentru forma a doua de apel). Aceste informatii vor fi exploatate de functiile ppval si unmkpp (a se vedea in continuare) .

Functia MATLAB ppval

Scop Calculeaza valorile aproximative ale unei functii aproximata polinomial pe subintervale.

Apel ye = ppval(pp,xe)

Parametri formali

Ø INTRARE : pp - matrice continand nodurile (intervalelor) pentru care s-a facut aproximarea, precum si coeficientii polinoamelor, pentru fiecare subinterval. Este obtinuta, de exemplu, cu functia splinen.

xe - valoare/ vector de valori pentru care se doreste aproximarea.

Ø IESIRE : ye - valorile aproximante calculate, corespunzatoare argumentelor din vectorul xe.

Functia MATLAB unmkpp

Scop Descompune o matrice de tip pp - de aproximare polinomiala pe portiuni - in componentele ei.

Apel [xi,c,l,k] = unmkpp(pp)

Parametri formali

Ø INTRARE : pp - matrice continand nodurile (intervalelor) pentru care s-a facut aproximarea, precum si coeficientii polinoamelor, pentru fiecare subinterval. Este obtinuta, de exemplu, cu functia splinen.

Ø IESIRE : xi - vector continand nodurile (intervalelor) pentru care s-a facut aproximarea.

c - matrice continand, pe fiecare linie, coeficientii polinoamelor pentru fiecare subinterval, in ordinea descrescatoare a puterilor.

l - numar de polinoame aproximante (= numar intervale = numar linii matrice c).

k - gradul polinoamelor de aproximare (= numar coloane matrice c).

Functia MATLAB mkpp

Scop Sinteza matrice de tip pp, pentru aproximarea polinomiala pe portiuni.

Apel pp = mkpp(xi,c)

Parametri formali

Ø INTRARE : xi - vector continand punctele (argumentele) pentru care s-a realizat aproximarea (interpolarea); are cel putin doua elemente.

c - matrice continand, pe fiecare linie, coeficientii polinoamelor, pentru fiecare subinterval definit de xi, in ordinea descrescatoare a puterilor. Are cel putin o linie. Conditie: numar linii c = numar elemente xi - 1 !

OBSERVATIE:

Pentru aproximarea cu functii spline cubice procedura de urmat este urmatoarea:

- aproximare valori functie:

ye = splinen(xi,yi,xe);

- aproximare valori functie, valori derivate, integrala definita:

pp = splinen(xi,yi);

ye = ppval(pp,xe); - aproximare functie

[xi,c,l,g] = unmkpp(pp);

c Þ cd1 (ecuatia (2)) - derivata I-a

ppd1 = mkpp(xi,cd1);

d1e = ppval(ppd1,xe);

c Þ cd2 (ecuatia (3)) - derivata a-II-a

ppd2 = mkpp(xi,cd2);

d2e = ppval(ppd2,xe);

c Þ cq (ecuatia (4)) - integrala definita pe intervalul [aq,bq];

aq=X(i); bq=X(i+1)

dx = bq - aq

q = polyval(cq,dx) - polyval explicata in continuare

Functia MATLAB polyval

Scop Evaluarea valorii/ valorilor unui polinom

Apel ye = polyval(c,xe)

Parametri formali

Ø INTRARE : c - vector continand coeficientii polinomului, in ordinea descrescatoare a puterilor.

xe - valoare/ vector de valori pentru care se doreste evaluarea polinomului.

Ø IESIREA : ye - valorile polinomului pentru argumentele furnizate in vectorul xe.

APLICATIE:

Problema

Se considera functia de test . Functia se evalueaza pe intervalul dat, cu un pas pentru interpolare . Se genereaza vectorii pentru interpolare xi, yi. Se determina functia spline cubica naturala corespunzatoare. Pentru testarea calitatii aproximarii, se evalueaza functia spline astfel obtinuta utilizand un pas de evaluare , he < hi, astfel generand vectorul xe. Se determina, apoi, vectorul ye ce contine valorile aproximante. De asemenea, se determina functiile spline pentru aproximarea derivatei de ordin I si derivatei de ordin al II-lea si se aproximeaza integrala definita pe un interval [aq,bq]. Se obtin astfel, respectiv, vectorii d1e, d2e si respectiv valoarea q. Valoarea q se compara cu cea teoretica:

Perechile de date: (xi,yi), (xe,ye), (xe,d1e), (xe,d2e), precum si: ([aq aq], [-1.25 1.25]), ([bq bq], [-1.25 1.25]) ce marcheaza intervalul pe care s-a facut integrarea, se vor reprezenta grafic. Imaginea grafica va avea inscrisa si valoarea aproximativa, q, a integralei.

Programul

Sa se realizeze, in mediul de programare MATLAB, programul pentru rezolvarea problemei expuse anterior, pentru testarea calitatii aproximarii numerice realizate cu ajutorul functiilor spline cubice naturale.

In realizarea programului se vor parcurge urmatoarele etape:

1. Introducere pas date de interpolat, in variabila hi; conditie de indeplinit: .

2. Generare vectori pentru interpolare, cu secventa de instructiuni MATLAB:

xi = 0:hi:4∙pi;

yi = sin(xi);

3. Generare matrice pp apeland functia splinen; secventa MATLAB:

pp = splinen(xi,yi)

Extragere informatii pentru derivate si integrala definita:

[x,c,l,g] = unmkpp(pp);

AI = c(:,1); BI = c(:,2);

CI = c(:,3); DI = c(:,4);

4. Introducere pas de evaluare in variabila he; conditie de indeplinit:

i si he < hi .

5. Generare vector pentru evaluare, xe, utilizand pasul de discretizare he. Se va proceda intr-o maniera similara celei in care s-a generat vectorul xi.

6. Evaluare functie spline cubica pentru valorile din vectorul xe. Valorile rezultate in vectorul ye:

ye = ppval(pp,xe);

7. Generare si evaluare functie spline pentru derivata I-a, cu secventa de instructiuni MATLAB:

cd1 = [3*AI 2*BI CI];

ppd1 = mkpp(xi,cd1);

d1e = ppval(ppd1,xe);

8. Generare si evaluare functie spline pentru derivata a II-a, cu secventa de instructiuni MATLAB:

cd2 = [6*AI 2*BI];

ppd2 = mkpp(xi,cd2);

d2e = ppval(ppd2,xe);

9. Introducere interval de integrare [aq,bq]; conditie de indeplinit:

.

10. Calcul integrala definita. Pentru inceput se determina intervalele unde se gasesc, respectiv, aq (indexul ia) si bq (indexul ib). Secventa care realizeaza acest lucru este descrisa in pseudocod:

ATRIBUIE k ← 0

ATRIBUIE P ← 1

_CAT TIMP ( P > 0 ) EXECUTA

| ATRIBUIE k ← k+1

|_ ATRIBUIE P ← (aq - xi(k)) ∙ (aq - xi(k+1))

ATRIBUIE ia ← k

ATRIBUIE k ← 0

ATRIBUIE P ← 1

_CAT TIMP ( P > 0 ) EXECUTA

| ATRIBUIE k ← k+1

|_ ATRIBUIE P ← (bq - xi(k)) ∙ (bq - xi(k+1))

ATRIBUIE ib ← k

ATRIBUIE dx ← aq - xi(ia)

_DACA ( dx ≠ 0 ) ATUNCI

| ATRIBUIE cq ← [(1/4)∙AI(ia) (1/3)∙BI(ia) (1/2)∙CI(ia) DI(ia) 0]

| ATRIBUIE qma ← polyval(cq,dx)

|ALTFEL

|_ ATRIBUIE qma ← 0

ATRIBUIE nq ← ib-ia

ATRIBUIE ind ← ia-1

ATRIBUIE qs ← 0

_PENTRU iq=1:nq EXECUTA

| ATRIBUIE ind ← ind+1

| ATRIBUIE dx ← xi(ind+1) - xi(ind)

| ATRIBUIE cq ←[(1/4)∙AI(ind) (1/3)∙BI(ind) (1/2)∙CI(ind) DI(ind) 0]

|_ ATRIBUIE qs ← qs + polyval(cq,dx)

ATRIBUIE dx ← bq - xi(ib)

_DACA ( dx ≠ 0 ) ATUNCI

| ATRIBUIE cq ← [(1/4)∙AI(ib) (1/3)∙BI(ib) (1/2)∙CI(ib) DI(ib) 0]

| ATRIBUIE qpb ← polyval(cq,dx)

|ALTFEL

|_ ATRIBUIE qpb ← 0

ATRIBUIE q ← qs + qpb - qma

ATRIBUIE qc ← cos(aq) - cos(bq)

SCRIE 'q =', q, 'qc =', qc

11. Reprezentare grafica: se repezinta in acelasi ecran grafic perechile de vectori:

(xi,yi) (caracter 'o')

(xe,ye) (linie continua '-')

(xe,d1e) (linie intrerupta '--')

(xe,d2e) (linie punctata ':')

([aq aq],[-1.25 1.25]) (linie-punct '_.')

([bq bq],[-1.25 1.25]) (linie-punct '_.')

Informatii referitoare la grafic - instructiuni MATLAB:

xlabel('xi, xe')

title(['yi:o ye:- d1:-- d2: ' num2str(q)])

12. Reluare vizualizare grafic, daca se doreste (punctul 11).

13. Reluare calcul integrala, pentru alt interval [aq,bq] daca se doreste (punctul 9).

14. Reluare cu alt pas de evaluare he, daca se doreste (punctul 4).

15. Reluare program cu alt pas de interpolare hi, daca se doreste (punctul 1).

OBSERVATII :

1. Descrierea anterioara a programului foloseste urmatoarele functii MATLAB: splinen, unmkpp, ppval, mkpp, polyval, sin, cos.

2. Programul se va rula pentru urmatoarele date de intrare:

I.) hi = 1; he = 0.1; aq = 1.57; bq = 3.14;

aq = 0.0; bq = 3.14;

II.) hi = 0.25; he = 0.025; aq = 1.57; bq = 3.14

aq=0.0; bq=3.14