Ajuste de Curvas em MATLAB

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

Ajuste de Curvas/Ajuste de Curvas marli.ppt
*
*
*
Ajuste de Curvas
Método dos Mínimos Quadrados
Profa. Marli
Disciplina ST462
*
*
*
Método dos mínimos quadrados
Neste caso, não se trata de interpolação, mas de ajustar uma função previamente escolhida a um conjunto de pontos.
Natureza do problema físico.
O método consiste em minimizar a soma dos quadrados das diferenças:
*
*
*
Introdução
Em muitos situações, conhece-se uma tabela de pontos (xi,yi), onde cada yi é obtido experimentalmente. 
E deseja obter expressões analíticas de uma dada curva y=f(x) que melhor se ajusta a esse conjunto de pontos.
*
*
*
Exemplo
Sabe-se que o número y de bactérias, por unidade de volume, existente em uma cultura após um determinado número x de horas, cresce exponencialmente com o aumento de x. 
Nesse caso número de bactérias cresce com o decorrer das horas na forma
O Problema consiste então, em determinar os valores mais aproximados dos parâmetros  e  desta exponencial. 
*
*
*
*
*
*
Ajuste a uma Reta
Mostraremos inicialmente, como ajustar um conjunto de dados a uma reta onde a e b são parâmetros a serem determinados. 
Nesse caso estamos interessados em minimizar a distancia de cada ponto (xi,yi), da tabela à cada ponto da (xi,a+bxi) da reta conforme figura.
*
*
*
Distancia de um ponto (xi,yi) à reta
*
*
*
 A distancia entre esses pontos é: 
E a soma dos quadrados dessas distancias: 
	
 sendo n o número de pontos
As candidatas a ponto de mínimo da função acima são aqueles para os quais são nulos as derivadas de q em relação a cada um dos parâmetros a e b respectivamente.
*
*
*
Tendo em vista que:
E que:
*
*
*
Obtemos o seguinte sistema de equação, denominado “equações normais” do problema cujas incógnitas são os parâmetros a e b da equação y = a + bx:
 Resolvendo o sistema obtemos os coeficientes a e b da reta:
*
*
*
*
*
*
Para ajustar por um polinômio de grau m, obteremos um sistema de (m+1) equações normais.
*
*
*
Exemplo: ajuste de curva à uma parábola
			p(x)=b0+b1x+b2x2
*
*
*
Tabela de cálculos
Plan1
			x			y			x2			x3			x4			xy			x2y
			1.0900			1.3500			1.1881			1.2950			1.4116			1.4715			1.6039
			1.2800			1.5800			1.6384			2.0972			2.6844			2.0224			2.5887
			1.3600			1.6800			1.8496			2.5155			3.4210			2.2848			3.1073
			1.4400			1.8500			2.0736			2.9860			4.2998			2.6640			3.8362
			1.6000			2.2300			2.5600			4.0960			6.5536			3.5680			5.7088
			1.6500			2.3800			2.7225			4.4921			7.4120			3.9270			6.4795
			8.4200			11.0700			12.0322			17.4817			25.7824			15.9377			23.3244
Plan2
			
Plan3
			
*
*
*
Montar e resolver o sistema linear
retirando os cálculos somatórios (S) da tabela
*
*
*
Muitos casos de funções matemáticas que representam de fenômenos físicos podem ser reduzidos ao caso linear.
Exemplos de algumas funções:
*
*
*
*
*
*
*
*
*
Função a ser ajustada: y = a + bx.
Deve-se obter o mínimo de através de:
*
*
*
Medida da qualidade do ajuste
A qualidade de um ajuste linear pode ser verificada em função do coeficiente de determinação r2 dado por:
*
*
*
Quanto mais próximo de um r2 estiver, melhor o ajuste.
r2 = 1 ajuste perfeito
0,8 < r2 < 1 ajuste bom
0,5 < r2 < 0,8 ajuste ruim
r2 < 0,5 sem correlação entre os xi e yi
*
*
*
Ajuste de Curvas/AjusteLinear.m
function [a b] = AjusteLinear(n,x,y) 
%Minimos Quadrados - ajuste de curva Linear
% y = a+bx 
%Dados de entrada n( numero de pontos (x,y)), x e y
% para isso temos que formar um sistema da forma 
% na + bSx = Sy
% aSx + bSx^2 = Sxy
% onde Sx: Somatorio dos valores de x(i), i=1,...,n 
% Sy: Somatorio dos valores de y(i), i=1,...,n
% Sx^2^: Somatorio dos valores de x(i)^2, i=1,...,n (Sx2)
% Sxy: Somatorio dos valores de x(i)y(i), i=1,...,n 
% Resolvendo o sistema temos que:
% a = (SySx^2 - SxSxy)/(nSx^2-(Sx)^2)
% b = (nSxy-SxSy)/(nSx^2-(Sx)^2)
x2= x.^2; % x2(i)= x(i)^2, i = 1,...,n
xy= x.*y; % xy(i) = x(i)y(i), i = 1,...,n
% Somatorias 
Sx=0; Sy=0; Sxy = 0; Sx2=0; 
for i=1:n
 Sx = Sx + x(i);
 Sy = Sy + y(i);
 Sxy = Sxy + xy(i);
 Sx2 = Sx2 + x2(i);
end
% determinando os coeficientes a e b
a = (Sy*Sx2-Sx*Sxy)/(n*Sx2-Sx^2);
b = (n*Sxy-Sx*Sy)/(n*Sx2-Sx^2);
disp('y = a + bx');[a b]
Ajuste de Curvas/AjusteLy_1abx.m
function [a b]= AjusteLy_1abx(n,x,y)
 % função y=1/(a+bx)
 % reducao aa linear
 %1/y = a + bx; v=1/y, aa=a, bb=b,w=x
 v=1./y;
 w=x;
 [a b] = AjusteLinear(n,w,v);
 disp('y = 1/(a+bx)');[a b]
Ajuste de Curvas/AjusteLy_abx.m
function [a b]= AjusteLy_abx(n,x,y)
 % função y=ab^x
 % reducao aa linear
 %lny = lna + x lnb; v=lny, aa=lna, bb=lnb,w=x
 v=log(y);
 w=x;
 [aa bb] = AjusteLinear(n,w,v);
 a = exp(aa); b=exp(bb);
 disp('y = ab^x');[a b]
Ajuste de Curvas/AjusteLy_aebx.m
function [a b]= AjusteLy_aebx(n,x,y)
 % função y=ae^bx
 % reducao aa linear
 %lny = lna + bx; v=lny, aa=lna, bb=b,w=x
 v=log(y);
 w=x;
 [aa bb] = AjusteLinear(n,w,v);
 a = exp(aa); b=bb;
 disp('y = ae^bx');[a b]
Ajuste de Curvas/AjusteParabola.m
function [b0 b1 b2] = AjusteParabola(n,x,y) 
%Minimos Quadrados - ajuste de curva aa uma parábola p(x)
% p(x)= y = b0+b1x+b2x2
%Dados de entrada n(número de pontos (x,y)), x e y
%
% para isso temos que formar um sistema da forma 
%
%b = [b0 b1 b2];
% nb0 + b1Sx + b2Sx^2 = Sy
% b0Sx + b1Sx^2 + b2Sx^3 = Sxy
% b0Sx^2 + b1Sx^3 + b2Sx^4 = S(x^2)y
% onde Sx: Somatorio dos valores de x(i), i=1,...,n 
% Sy: Somatorio dos valores de y(i), i=1,...,n
% Sx^2^: Somatorio dos valores de x(i)^2, i=1,...,n (Sx2)
% Sx^3^: Somatorio dos valores de x(i)^3, i=1,...,n (Sx3) 
% Sx^4^: Somatorio dos valores de x(i)^4, i=1,...,n (Sx4)
% Sxy: Somatorio dos valores de x(i)y(i), i=1,...,n 
% S(x^2)y: Somatorio dos valores de (x(i)^2)y(i),i=1,...,n (Sx2y)
x2= x.^2; % x2(i)= x(i)^2, i = 1,...,n
x3= x.^3; % x3(i)= x(i)^3, i = 1,...,n
x4= x.^4; % x4(i)= x(i)^4, i = 1,...,n
xy= x.*y; % xy(i) = x(i)y(i), i = 1,...,n
x2y= (x.^2).*y; %x2y = (x(i)^2)y(i), i = 1,...,n
% Somatorias 
Sx=0; Sy=0; Sxy = 0; Sx2=0; Sx3=0; Sx4=0 ;Sx2y=0;
for i=1:n
 Sx = Sx + x(i);
 Sy = Sy + y(i);
 Sxy = Sxy + xy(i);
 Sx2 = Sx2 + x2(i);
 Sx3 = Sx3 + x3(i);
 Sx4 = Sx4 + x4(i);
 Sx2y = Sx2y + x2y(i);
end
% formar sistema Ab=d
% n Sx Sx2 b0 Sy
% A = Sx Sx2 Sx3 , b = b1 , d= Sxy
% Sx2 Sx3 Sx4 b2 Sx2y 
A = [n Sx Sx2;Sx Sx2 Sx3;Sx2 Sx3 Sx4];
d = [Sy; Sxy; Sx2y];
% Resolver o sistema b = (A^-1)d;
b= inv(A)*d;
disp('y = b0 + b1x + b2x^2')
b0=b(1); b1=b(2); b2=b(3);
[b0 b1 b2]
Ajuste de Curvas/CoefDeterm.m
function [ r2 ] = CoefDeterm(n,x,y,b0,b1,b2)
% Deterninar o coeficiente de determinacao r^2 do ajuste de curva 
% pelo métodos dos minimos quadrados
%
% Determinar a media de y medy= soma de y(i), i=1,...n 
medy = 0;
for i = 1:n 
 medy = medy + y(i);
end
medy = medy/n;
% r^2 = (soma de (p(x(i))- medy)^2, i=1,...,n)/(soma de (y(i)- medy)^2,
% i=1,...,n)
% determinar as Sp = (soma de (p(x(i))- medy)^2, i=1,...,n) e Sy =(soma de (y(i)- medy)^2,
% i=1,...,n)
Sp = 0; Sy =
0;
for i = 1:n
 Sp = Sp + (parabola(b0,b1,b2,x(i))-medy)^2;
 Sy = Sy + (y(i)-medy)^2;
end
r2 = Sp/Sy;
if r2 == 1
 disp('ajuste perfeito')
else
 if (r2 < 1) && (r2 >= 0.8)
 disp('ajuste bom')
 else
 if (r2 < 0.8) && (r2 >= 0.5)
 disp('ajuste ruim')
 else
 disp('sem correlacao entre x e y')
 end
 end
end
disp('coeficiente de determinação r2=');r2
 
 
Ajuste de Curvas/CoefDetermL.m
function [ r2 ] = CoefDetermL(n,x,y,yy)
% Deterninar o coeficiente de determinacao r^2 do ajuste de curva 
% pelo métodos dos minimos quadrados
%Dados de entrada 
% n: (numeros de pontos (x,y)), x , y 
% x e y: as coordenadas de (x,y)
% yy = f(x) (função o qual os pontos (x,y) foram ajustados) 
%
% Determinar a media de y medy= soma de y(i), i=1,...n 
medy = 0;
for i = 1:n 
 medy = medy + y(i);
end
medy = medy/n;
% r^2 = (soma de (p(x(i))- medy)^2, i=1,...,n)/(soma de (y(i)- medy)^2,
% i=1,...,n)
% determinar as Sp = (soma de (p(x(i))- medy)^2, i=1,...,n) e Sy =(soma de (y(i)- medy)^2,
% i=1,...,n)
Sp = 0; Sy = 0;
for i = 1:n
 Sp = Sp + (yy(i)-medy)^2;
 Sy = Sy + (y(i)-medy)^2;
end
r2 = Sp/Sy;
if r2 == 1
 disp('ajuste perfeito ')
else
 if (r2 < 1) && (r2 >= 0.8)
 disp('ajuste bom')
 else
 if (r2 < 0.8) && (r2 >= 0.5)
 disp('ajuste ruim')
 else
 disp('sem correlacao entre x e y')
 end
 end
end
disp('Coeficiente de Determinação r^2=');r2
 
 
Ajuste de Curvas/linear.m
function [ y ] = linear(a,b,x)
 % p(x) = b0 + b1*x + b2*(x^2);
 y = a + b.*x ;
Ajuste de Curvas/MinQuaLinear.m
clear all
clc
%Minimos Quadrados - ajuste de curva Linear
% y = a+bx 
%Dados de entrada n( numero de pontos (x,y)), x e y
n=6;
x = [1.0900 1.2800 1.3600 1.4400 1.6000 1.6500];
%x = [1.0900 1.0800 1.600 10.4400 1.6000 3.6500];
y = [1.3500 1.5800 1.6800 1.8500 2.2300 2.3800];
%Ajuste Linear
[a b]=AjusteLinear(n,x,y);
% coeficiente de determinacao
yy=linear(a,b,x);
r2= CoefDetermL(n,x',y',yy);
% Plotar Grafico
xx= min(x):.01:max(x);
yy = linear(a,b,xx);
%yy = a + b*xx;
plot(x,y,'r.',xx,yy,'b')
 
 
 
 
Ajuste de Curvas/MinQuaLinearRed.m
clear all
clc
%Minimos Quadrados - ajuste de curva Linear
% y = a+bx 
%Dados de entrada n( numero de pontos (x,y)), x e y
n=6;
x = [1.0900 1.2800 1.3600 1.4400 1.6000 1.6500];
y = [1.3500 1.5800 1.6800 1.8500 2.2300 2.3800];
%Ajuste de funcoes redutiveis aa Linear
[a b]= AjusteLy_abx(n,x,y); % Ajuste y=ab^x
% coeficiente de determinacao
yy=y_abx(a,b,x);
r2= CoefDetermL(n,x',y',yy);
% Plotar Grafico
xx= min(x):.01:max(x);
yy = y_abx(a,b,xx); %y=ab^x
plot(x,y,'r.',xx,yy,'b')
hold on 
%Ajuste de funcoes redutiveis aa Linear
[a b]= AjusteLy_aebx(n,x,y); % Ajuste y=ae^bx
% coeficiente de determinacao
yy=y_aebx(a,b,x);
r2= CoefDetermL(n,x',y',yy);
% Plotar Grafico
xx= min(x):.01:max(x);
yy = y_aebx(a,b,xx); %y=ae^bx
plot(xx,yy,'m')
hold on
%Ajuste de funcoes redutiveis aa Linear
[a b]= AjusteLy_1abx(n,x,y); % Ajuste y=1/(a+bx)
% coeficiente de determinacao
yy=y_1abx(a,b,x);
r2= CoefDetermL(n,x',y',yy);
% Plotar Grafico
xx= min(x):.01:max(x);
yy = y_1abx(a,b,xx); %y=1/(a+bx)
plot(xx,yy,'g')
 
 
 
 
Ajuste de Curvas/MinQuaParabola.m
clear all
clc
%Minimos Quadrados - ajuste de curva aa uma parábola p(x)
% p(x)= y = b0+b1x+b2x2
%Dados de entrada n(número de pontos (x,y)), x e y
n=6;
x = [1.0900 1.2800 1.3600 1.4400 1.6000 1.6500];
%x = [1.0900 1.0800 1.600 1.4400 1.6000 1.6500];
y = [1.3500 1.5800 1.6800 1.8500 2.2300 2.3800];
%
[b0 b1 b2] = AjusteParabola(n,x,y);
% coeficiente de determinação
r2=CoefDeterm(n,x',y',b0,b1,b2);
% Plotar Grafico
xx= min(x):.01:max(x);
yy = parabola(b0,b1,b2,xx);
%yy=b0 + b1*xx + b2.*(xx.^2);
plot(x,y,'r.',xx,yy,'b')
 
 
 
 
Ajuste de Curvas/parabola.m
function [ y ] = parabola(b0,b1,b2,x)
 % p(x) = b0 + b1*x + b2*(x^2);
 y = b0 + b1.*x + b2.*x.^2;
Ajuste de Curvas/y_1abx.m
function [y] = y_1abx(a,b,x)
 % y = 1/a+bx
 y = 1./(a+b.*x);
 
Ajuste de Curvas/y_abx.m
function [y] = y_abx(a,b,x)
 % y = ab^x
 y = a*b.^x;
 
Ajuste de Curvas/y_aebx.m
function [y] = y_aebx(a,b,x)
 % y = ae^bx
 y = a*exp(b.*x);