CURSUL 3. APROXIMAREA FUNCIILOR PRIN

AJUSTARE

Ajustarea datelor prin Metoda Celor Mai Mici Ptrate

Fie n+1 puncte x0, x1, .., xn n care se cunosc valorile unei funcii f(x). Metoda

celor mai mici ptrate se folosete la prelucrarea matematic a datelor x0, x1, .., xn

n urmtoarele cazuri: cnd n este foarte mare, sau cnd valorile funciei n

punctele considerate nu sunt exacte, fiind obinute experimental.

Prin aceast metod funcia f(x) se aproximeaz cu

,)(*)(

0====

====m

i

iim xcxF cu m0

s fie minim. nlocuim Fm(x) cu o funcie de constantele c0, .., cm. Se pune

problema determinrii acestor constante cj, j=0..m, astfel nct

==== ====

====n

i

m

j

ijjim xcfcc

0 0

20 )]([),..,(

s fie minim. Determinarea constantelor se face rezolvnd sistemul

mici

..0,0 ========

adic a sistemului

0)(])([ 00 0

==== ==== ====

i

n

i

m

j

ijji xxcf

0)(])([

0 0

==== ==== ====

im

n

i

m

j

ijji xxcf

Sistemul de mai sus se mai poate scrie i sub forma urmtoare:

==== ========

====n

i

n

i

ii

m

j

iijj xfxxc

0 00

00 )()()(

==== ========

====n

i

n

i

imi

m

j

imijj xfxxc

0 00

)()()(

De obicei n practic se consider mixx ii ..0,)( ======== , caz n care avem

Fm(x)=c0+c1x++cmxm,

iar sistemul de rezolvat este:

============

====++++++++++++++++n

i

i

n

i

mim

n

i

i fxcxcnc

00010 ..)1(

..

============

++++

====

====++++++++++++n

i

mii

n

i

mim

n

i

mi

n

i

mi xfxcxcxc

00

2

0

11

00 ..

Exemplu:

S considerm tabelul:

X 1 2 3 4 5 6 7

Y 2 2,5 3,5 6 6,5 8 9

Se vede c datele se aeaz aproximativ de-a lungul unei drepte y=ax+b. Dreapta,

dat de metoda celor mai mici ptrate, care ajusteaz mulimea dat de puncte are

coeficienii a i b dai de condiia ca expresia ( )E a b y ax bi ii

( , ) = 2 s fie

minim.

Aceasta implic:

E a b

a

E a b

b

( , )

( , )

=

=

0

0

(1)

Se obine sistemul:

S a S b S

S a pb S

x x xy

x y

2 1

1 1

+ =

+ =

(2)

unde :

S xx ii

22= S xx i

i1 = S x yxy i i

i

= S yy i1 = (3)

iar p este numrul de date. Avem mai jos graficul acestei drepte.

y = 1.25x + 0.3571

0

1

2

3

4

5

6

7

8

9

10

1 2 3 4 5 6 7

Mai general putem considera ajustari ale lui f prin expresii de forma

f(x)a0+a1x+a2x2+....+ahx

h. Coeficienii ai, conform metodei celor mai mici

ptrate sunt dai de condiia ca expresia E a a a y a a xh k hh

k p

( , , ... ) ( ... )..

0 1 02

0

= = s fie

minim, deci:

E

ai h

i

= =0 0 1, , ,... (4)

Rezolvarea sistemului (4) este dificila n multe situatii, sistemul fiind n general

neliniar si uneori singular. In cazul cnd sistemul este liniar rezolvarea lui este

relativ simpla.

Algoritmul Levenberg-Marquardt

Chiar daca sistemul (4) este singular functia E a a a y a a xh k hh

k p

( , , ... ) ( ... )..

0 1 02

0

= =

are un minim si programele standard de matematica precum scilab sau matlab l

determina. Sistemul (4) se mai poate scrie:

h,...,kpentru,NaMk

h

lll,k

211

===

unde + ==i

i

k

iki

lk

il,kyxNsixM 12 .

Sunt situatii cand modelul ales este neadecvat; Un program care sa determine o

formula potrivita pentru ajustarea unor date experimentale ar trebui sa aiba

memorate un numar mare de modele de functii model depinznd de unul sau mai

multi parametri si o metoda de a determina minimul expresiei

E a a a y a a xh k hh

k p

( , ,... ) ( ... )..

0 1 02

0

= = pentru fiecare model.

Modelul pentru care =

=n

iini)y)a,...,a,x(f(

nR

1

2

1

1 are cea mai mica valoare este

n principiu modelul cel mai bun. Divizarea cu n este facuta pentru a anula efectul

numarului de puncte asupra lui R.

Pentru functii f : Rm R depinznd de ( )m

x,...,x,xx 21=r

tehnica ajustarii este

aceeasi, doar ca functiile model depind de mRxr

, deci sunt de forma

)a,...,a,x(fp1

r. Avem

( )211 =i

ipipy)a,...,ax(f)a,...,a(E

r (5)

unde mm,i,ii

R)a,...,x(x = 1r

Pentru minimizarea lui )a,...a,a(Ep21 exista un algoritm numit "Levenberg-

Marquardt" care este implementat n mai multe programe ca, de exemplu, scilab

sau matlab.

Minimalizarea expresiei (5) in scilab

Ca date de intrare avem un tabel cu date experimentale n..iii

)y,x( 1=r

si o functie

model )a,...,a,x(fp1

r. Se cere determinarea parametrilor

pa,...,a1 astfel ca (5) sa

fie minima. Procedura scilab pentru determinarea parametrilor este lsqrsolve care are forma

[a,v,[info]]=lsqrsolve(a0,e,n,[stop],[diag])

unde:

a - vectorul cu componentele cautate p

a,...,a1 .

v - vectorul cu diferentele ipi

y)a,...,ax(f 1r

pentru valoarea lui a gasita de program

info - este numar a carui valoare depinde de modul de terminare al programului: 0

- parametrii de intrare sunt incorecti; 1-algoritmul estimeaza atol (vezi intrarea

stop), adica eroarea relativa ntre a si solutia exacta, a fost atinsa; 2 - a fost atins

numarul maxim maxeval (vezi intrarea stop) de evaluari pentru functia e; 3 - atol

este o valoare prea mica si nu poate fi atinsa; 4 - iteratiile nu progreseaza satisfacator.

a0 - un vector cu valori initiale pentru a

e - o functie de forma e(a, n) unde a este vectorul parametrilor p

a,...,a1 , iar n este

numarul de date ( )ii

y,xr

din tabelul de intrare. Valoarea functiei e este un vector cu

n componente

( )( )

( )

nnya,xf

..................

ya,xf

ya,xf

r

r

r

22

11

stop - este un vector optional cu conduitii de oprire a algoritmului, stop=[etol,

atol, gol, maxeval, epsfcn, factor], cu valorile implicite [1.d-8,1.d-8,1.d-

5,1000,0,100].

Semnificatiile unor parametri sunt:

atol - eroarea relativa maxima admisa intre doua valori iterative consecutive ale lui

a (daca se obtine o eroare relativa mai mica, programul se opreste)

maxeval - numarul maxim de apeluri ale functiei e (daca se atinge acest maxim,

programul se opreste)

diag - este un vector de ponderi

Un exemplu de utilizare a acestei rutine este dat mai jos.

a1=1;a2=-1;a3=2;a4=0.5; a0=[10;10;-10;10]; a=a0; //definim o functie model deff(y=F(x),y=a1*x(1)+a2*x(2)+a3*sin(x(3))+a4*x(4)*x(5)); //definim o matrice de date data=rand(500,5); //definim y=F(x) y=zeros(500,1); for i=1:500 x=data(i,1:5); y(i)=F(x); end //stop=[1.d-10, 1.d-12, 1.d-5, 50000, 0, 100]; //ne pregatim de optimizare //functia de optimizat function dif=e(a,n) dif=a(1)*data(:,1)+a(2)*data(:,2)+a(3)*sin(data(:,3))... +a(4)*data(:,4).*data(:,5)-y; endfunction //este apelataa rutina de optimizare [a,v]=lsqrsolve(a0,e,size(data,1)); a norm(v,.inf.)

Rezultatele ntoarse sunt a = ! 1. ! ! - 1. ! ! 2. ! ! 0.5 !

Valorile exacte ale parametrilor sunt 1, -1, 2, 0.5. Norma vectorului diferentelor

ipiy)a,...,ax(f 1

r este 4.441D-16 .

In programul de mai sus am construit pe y dupa formula

544332211 ,i,i,i,i,ii xxa)xsin(axaxay +++= si am cerut programului sa determine (a1, a2, a3, a4) plecnd de la a0 = (10, 10;-10; 10) n functia model ( ) 5443322114321 xxa)xsin(axaxaa,a,a,a,xf +++=r

. Dupa cum se vede s-a obtinut valoarea exacta pentru a.

Pentru a apela procedura pentru un numar mai mare de modele putem aranja calculele ca mai jos:

// multitest pentru functii de o variabila //valoarea initiala a parametrilor a a0=[1;1;1;1;1;1;1;1;1;1;1]; //numarul maxim admis de parametri necunoscuti este 10; //functiiile model pentru ajustare nrmodele=3; function val=model(itip,a,x) select itip case 1 then val=a(1)+a(2)*x; case 2 then val=a(1)+a(2)*x+a(3)*x.^2; case 3 then val=a(1)+a(2)*x+a(3)*x.^2+a(3)*x.^3; end endfunction //functia pentru procedura lsqrsolve function dif=e(a,n); dif=model(itip,a,x)-y; endfunction //datele de ajustat x=[1;2;3;4;5;6;7;8;9;1.5;2.5;3.5]; y=[2;5;9;17;24;40;50;65;80;3;5;10]; n=length(x); // testarea tuturor modelelor r=zeros(nrmodele,1); for itip=1:nrmodele [a,v]=lsqrsolve(a0,e,n); //discrepanta intre model si datele reale r(itip)=norm(v,2)/n; end //discrepanta minima [rmin,itip]=min(r); //parametrii la discrepanta minima [a,v]=lsqrsolve(a0,e,n); //tiparim punctele si functia care ajusteaza plot(x,y,ro); xgrid; stanga=min(x);dreapta=max(x); t=(stanga:0.1:dreapta); vt=model(itip,a,t); plot(t,vt,b:);

Se obtine urmatorul grafic, pe care se vad punctele date si graficul functiei care le ajusteaza. Dintre cele trei modele, n acest caz modelul 2 este cel mai bun (r este minim)

Exemplu

S se ajusteze printr-un polinom de gradul m=3, urmtoarele date experimentale:

x=[-2, -1, 0, 1, 2] i y=[12, 6, 2, 0, 0]

Rezolvare n Scilab

Program de ajustare printr-un polinom generalizat de grad m. Apelarea se face

pentru datele din exemplu.

function c=CeleMaiMiciPatrate(x,y,m) n=length(x); for k=1:2*m+1 s(k)=0; for i=1:n s(k)=s(k)+x(i)^(k-1); end end for k=1:m+1 b(k)=0; for i=1:n b(k)=b(k)+y(i)*x(i)^(k-1); end end for i=1:m+1 for j=1:m+1 a(i,j)=s(i+j-1); end

end c=a^(-1)*b; endfunction CeleMaiMiciPatrate([-2 -1 0 1 2],[12 6 2 0 0],3) ans = ! 2. ! ! - 3. ! ! 1. ! ! 0. !

Exerciii

1. S se scrie expresia polinomului din exemplul de mai sus i s se reprezinte pe

acelai grafic datele experimentale i polinomul.

2. Se consider urmtoarele date experimentale:

xi 0 0.5 1 1.5

yi 1 1.1276 1.5431 2.3534

Se cere:

a) s se calculeze valoarea lui y pentru x=0.25 folosind o parabol pentru ajustarea

datelor experimentale prin metoda celor mai mici ptrate;

b) s se ajusteze aceste date printr-o hiperbol folosind metoda celor mai mici

ptrate, apoi s se reprezinte grafic;

c) s se ajusteze prin y=a*sin(x)+b datele experimentale de mai sus.