UNIVERSIDAD AUTÓNOMA DE YUCATÁN
FACULTAD DE INGENIERÍA QUÍMICA
MÉTODOS NUMÉRICOS
ADA 6: polinomios de interpolación y ajuste de curvas
Elaborado por:
Cetz Guzmán Ángel Ramiro
Cob Pérez Marco Antonio
Mis Navarrete Isaac Uriel
Solís Montero Nehemías Aarón
Tabasco Matos Argel Antonio
Fecha de entrega:
15 de noviembre, 2019
Ejercicio 1. Dados los datos
𝑥 1 2 2.5 3 4 5
𝑓(𝑥) 1 5 7 8 2 1
estime el valor de 𝑓(3.4) usando interpolación lineal y luego con interpolación
cuadrática.
Interpolación lineal
P1=input('ingrese primer par de datos[a,f(a)]\n');
P2=input('ingrese primer par de datos[b,f(b)]\n');
x=input('Ingrese el punto a aproximar x\n');
a= P1(1);
fa=P1(2);
b=P2(1);
fb=P2(2);
y=(fb-fa)/(b-a)*(x-a)+fa;
fprintf('f(%f)=%f',x,y);
>> x=[1 2 2.5 3 4 5];
>> y=[1 5 7 8 3 2];
>> inter_lineal
ingrese primer par de datos[a,f(a)]
[3,8]
ingrese primer par de datos[b,f(b)]
[4,2]
Ingrese el punto a aproximar x
3.4
f(3.400000)=5.600000
Interpolación cuadrática
%Entradas
P1=input('ingrese primer par de datos[a,f(a)]\n');
P2=input('ingrese segundo par de datos[b,f(b)]\n');
P3=input('ingrese tercer par de datos[c,f(c)]\n');
x=input('Ingrese el punto a aproximar x\n');
%Asignacion de valores
a= P1(1);
fa=P1(2);
b=P2(1);
fb=P2(2);
c=P3(1);
fc=P3(2);
% Calculo de coeficiente
a0=fa;
a1=(fb-fa)/(b-a);
a2=((fc-fb)/(c-b)-(fb-fa)/(b-a))/(c-a);
%Salida
y=a0+a1*(x-a)+a2*(x-a)*(x-b);
fprintf('f(%f)=%f\n',x,y)
>> inter_cuad
ingrese primer par de datos[a,f(a)]
[3,8]
ingrese segundo par de datos[b,f(b)]
[4,2]
ingrese tercer par de datos[c,f(c)]
[5,1]
Ingrese el punto a aproximar x
3.4
f(3.400000)=5.000000
Ejercicio 2. Dados los datos
𝑥 1 2 2.5 3 4 5
𝑓(𝑥) 1 5 7 8 2 1
1) Calcule 𝑓(3.4) usando el polinomio interpolador de Newton
2) Estime el error de su predicción.
>> x=[1 2 2.5 3 4 5];
>> y=[1 5 7 8 2 1];
function [D,yi,p,xi]=internw(x,y,xi)
%Input - x es un vector que contiene una lista de abcisas
% - y es un vector que contiene una lista de ordenadas
% - xi es el valor a hallar
%Output - yi es el valor hallado
% - D es la tabla de diferencias de diviciones
n=length(x);
D=zeros(n,n);
D(:,1)=y';
%Formula para la tabla de diferencia de divisiones
for j=2:n
for k=j:n
D(k,j)=(D(k,j-1)-D(k-1,j-1))/(x(k)-x(k-j+1));
end
end
%Inicializa las variables
n=length(x);
b=zeros(n);
b(:,1)=y(:);
%Obtener la tabla de diferencias
for j=2:n
for i=1:n-j+1
b(i,j)=(b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
end
end
%Calcular el dato interpolado
xt=1;
yi=b(1,1);
for j=1:n-1
xt=xt.*(xi-x(j));
yi=yi+b(1,j+i)*xt;
end
%Construir el polinomio
p=num2str(b(1,1));
xx=x*-1;
for j=2:n
signo='';
if b(1,j)>=0
signo='+';
end
xt='';
for i=1:j-1
signo2='';
if xx(i)>=0
signo2='+';
end
xt=strcat(xt,'*(x',signo2,num2str(xx(1)),')');
end
p=strcat(p,signo,num2str(b(1,j)),xt);
end
%'para calcular el error es D*(xi-x1)*(xi-x2)*(xi-xn)'
end
>> xi=3.4
xi =
3.4000
>> [D,yi,p,xi]=internw(x,y,xi)
D=
Columns 1 through 5
1.0000 0 0 0 0
5.0000 4.0000 0 0 0
7.0000 4.0000 0 0 0
8.0000 2.0000 -2.0000 -1.0000 0
2.0000 -6.0000 -5.3333 -1.6667 -0.2222
1.0000 -1.0000 2.5000 3.1333 1.6000
Column 6
0.4556
yi =
6.9766
p=
1+4*(x-1)+0*(x-1)*(x-1)-1*(x-1)*(x-1)*(x-1)-0.22222*(x-1)*(x-1)*(x-
1)*(x-1)+0.45556*(x-1)*(x-1)*(x-1)*(x-1)*(x-1)
xi =
3.4000
Para calcular el error:
>>Error= 0.4556*(3.4-1)*(3.4-2)*(3.4-2.5)*(3.4-3)*(3.4-4)*(3.4-5)
Error =
0.5291
Ejercicio 3. Dados los datos
𝑥 1 2 3 5 6
𝑓(𝑥) 4.75 4 5.25 19.75 36
Calcule 𝑓(4) usando el polinomio interpolador de Lagrange
function [C,L]= lagran(X,Y)
W=length(X);
n= W-1;
L= zeros(W,W);
%Formando los coeficientes de lagrange
for k = 1: n+1
V= 1;
for j =1: (n+1)
if k~= j
V= conv(V, poly(X(j)))/(X(k)-X(j));
end
end
L(k,:)=V;
end
%Coeficientes del polinomio
C=Y*L;
X=[1 2 3 5 6];
>> Y=[4.75 4 5.25 19.75 36];
>> [C,L]= lagran(X,Y)
C=
0.0000 0.2500 -0.5000 -1.0000 6.0000
L=
0.0250 -0.4000 2.2750 -5.4000 4.5000
-0.0833 1.2500 -6.4167 12.7500 -7.5000
0.0833 -1.1667 5.4167 -9.3333 5.0000
-0.0417 0.5000 -1.9583 3.0000 -1.5000
0.0167 -0.1833 0.6833 -1.0167 0.5000
>> polyval(C,4)
ans =
10.0000
Ejercicio 4. En estudios de polimerización inducida por radiación, se emplea
una fuente de rayos gamma para obtener dosis medidas de radiación. Sin
embargo, la dosis varía con la composición del aparato según los datos que se
muestran a continuación:
Posición
en pulgadas Dosis en 105 rads/h
desde el punto base
0 1.9
0.5 2.39
1 2.71
1.5 2.98
2 3.2
3 3.2
3.5 2.98
4 2.74
Como puede observarse, la lectura en 2.5 pulgadas no se reportó, pero se
necesita el valor de esta radiación. Utilice polinomios de interpolación de varios
grados a los datos para estimar el dato faltante. ¿Algún polinomio de mínimos
cuadrados proporciona mejor ajuste? Justifique su respuesta.
>> x=[0 0.5 1 1.5 2 3 3.5 4];
>> y=[1.9 2.39 2.71 2.98 3.2 3.2 2.98 2.74];
>> C=polyfit(x,y,7)
Columns 1 through 5
-0.0024 0.0312 -0.1385 0.2115 0.0858
Columns 6 through 8
-0.6348 1.2572 1.9000
La ecuación encontrada es y= -0.0024x7+0.0312x6-
0.1385x5+0.2115x4+0.0858x3-0.6348x2+1.2572x+1.9000
>> polyval(C,2.5)
ans =
3.2907
El valor 3.2907 es la dosis en 105rads/h cuando la posición en pulgadas desde
el punto base es de 2.5
Mínimos cuadrados
Usaremos las siguientes fórmulas para calcular a0 y a1
>> Exi=sum(x)
Exi =
15.5000
>> Eyi=sum(y)
Eyi =
22.1000
>> Exiyi=sum(x.*y)
Exiyi =
45.7650
>> Exi2=sum(x.^2)
Exi2 =
44.7500
>> a1=(8*Exiyi-(Eyi*Exi))/((8*Exi2)-(Exi.^2))
a1 =
0.2002
>> a0=(Eyi-(a1*Exi))/8
a0 =
2.3747
Ecuación de la recta es y= 2.3747+0.2002x
>> x1=linspace(0,4);
>> y1=2.3747+0.2002*x1;
>> plot(x,y,'o',x1,y1,'r'),grid
Polinomio de 2do grado
Usaremos el siguiente método para calcular el polinomio de grado 2
Conocemos algunos datos, pero se debe calcular los datos faltantes para poder
encontrar los valores de a0, a1 y a2
>> Ex3=sum(x.^3)
Ex3 =
146.3750
>> Ex4=sum(x.^4)
Ex4 =
509.1875
>> Exi2yi=sum(y.*x.^2)
Exi2yi =
131.9575
Utilizando el método de la matriz inversa para encontrar los valores
>> A=[8 Exi Exi2; Exi Exi2 Ex3; Exi2 Ex3 Ex4]
A=
8.0000 15.5000 44.7500
15.5000 44.7500 146.3750
44.7500 146.3750 509.1875
>> A=[8 Exi Exi2; Exi Exi2 Ex3; Exi2 Ex3 Ex4]
A=
8.0000 15.5000 44.7500
15.5000 44.7500 146.3750
44.7500 146.3750 509.1875
>> B=[22.10;45.7650;131.9575]
B=
22.1000
45.7650
131.9575
>> X=A\B
X=
1.8932
1.0634
-0.2129
Lo que significa que los valores son:
a0=1.8932
a1=1.0634
a2=-0.2129
La ecuación es y=1.8932+1.0634x-0.2129x2
>> x1=linspace(0,4);
>> y1=1.8932+1.0634*x1-0.2129*x1.^2;
>> plot(x,y,'o',x1,y1,'r'),grid
Polinomio de 3er grado
Calculamos los datos faltantes
>> Ex5=sum(x.^5)
Ex5 =
1.8328e+03
Ex6 =
6.7397e+03
Exi3yi =
428.1938
Utilizando el método de la matriz inversa para encontrar los valores
>> A=[8 Exi Exi2 Ex3; Exi Exi2 Ex3 Ex4; Exi2 Ex3 Ex4 Ex5;Ex3 Ex4 Ex5 Ex6]
A=
1.0e+03 *
0.0080 0.0155 0.0447 0.1464
0.0155 0.0447 0.1464 0.5092
0.0447 0.1464 0.5092 1.8328
0.1464 0.5092 1.8328 6.7397
>> B=[22.10;45.7650;131.9575;428.1938]
B=
22.1000
45.7650
131.9575
428.1938
>> X=A\B
X=
1.9061
1.0054
-0.1727
-0.0068
Lo que significa que los valores son:
a0=1.9061
a1=1.0054
a2= -0.1727
a3= -0.0068
La ecuación es y=1.9061+1.0054x-0.1727x2-0.0068x3
>> x1=linspace(0,4);
>> y1=1.9061+1.0054*x1-0.1727*x1.^2-0.0068*x1.^3;
>> plot(x,y,'o',x1,y1,'r'),grid
Evaluando las funciones para 2.5 se obtiene:
>> x1=2.5;
Polinomio de grado 2
>> y1=1.8932+1.0634*x1-0.2129*x1.^2
y1 =
3.2211
Polinomio de grado 3
>> y1=1.9061+1.0054*x1-0.1727*x1.^2-0.0068*x1.^3
y1 =
3.2340
Como se puede observar en las gráficas y en las evaluaciones, el ajuste se va
aproximando al valor original que se obtuvo con mínimos cuadrados
Ejercicio 5. Considere el siguiente conjunto de datos
𝑥 1 2 3 4 5
𝑓(𝑥) 0.6 1.9 4.3 7.6 12.6
determine la curva que mejor se ajusta en el sentido de mínimos cuadrados,
llevando a cabo los siguientes pasos:
1) Ajuste una curva de la forma 𝑓(𝑥) = 𝐶𝑒 𝐴𝑥
2) Ajuste una curva de la forma 𝑓(𝑥) = 𝐶𝑥 𝐴
3) Use el error cuadrático medio dado por:
1
𝑛 2
1
𝐸2 (𝑓) = ( ∑|𝑓(𝑥𝑘 ) − 𝑦𝑘 |2 )
𝑛
𝑘=1
para determinar cuál de las dos curvas ajusta mejor.
>> x=[1 2 3 4 5];
>> y=[0.6 1.9 4.3 7.6 12.6];
>> Exi=sum(x)
Exi =
15
>> Eyi=sum(y)
Eyi =
27
>> Exiyi=sum(x.*y)
Exiyi =
110.7000
>> Exi2=sum(x.^2)
Exi2 =
55
>> a1=(5*Exiyi-(Eyi*Exi))/((5*Exi2)-(Exi.^2))
a1 =
2.9700
>> a0=(Eyi-(a1*Exi))/5
a0 =
-3.5100
Ecuación de la recta es y= -3.5100+2.9700x
1) Ajuste una curva de la forma 𝑓(𝑥) = 𝐶𝑒 𝐴𝑥
>> x=[1 2 3 4 5];
>> y=[0.6 1.9 4.3 7.6 12.6];
>> Y=log(y)
Y=
-0.5108 0.6419 1.4586 2.0281 2.5337
>> polyfit(x,Y,1)
ans =
0.7475 -1.0123
La recta es: y= 0.7475x-1.0123
Modelo exponencial es:
>> y=exp(-1.0123)*exp(0.7475*x);
>> X=linspace(0,5);
>> Y=-3.5100+2.9700*X;
>> x1=linspace(0,5);
>> y1=exp(-1.0123)*exp(.7475*x1);
>> plot(X,Y,'b',x1,y1,'g--'),grid
>> legend('Funcion Minimos cuadrados','y=Ce^Ax')
Error
>> x=[1 2 3 4 5];
>> y=-3.5100+2.9700*x;
>> y1=exp(-1.0123)*exp(.7475*x);
>> Er=(1/5*((sum(y-y1)).^2)).^(1/2)
Er =
0.5793
2) Ajuste una curva de la forma 𝑓(𝑥) = 𝐶𝑥 𝐴
>> x=[ 1 2 3 4 5];
>> y=[0.6 1.9 4.3 7.6 12.6];
>> X=log(x);
>> Y=log(y);
>> polyfit(X,Y,1)
ans =
1.8860 -0.5755
La recta es: 1.8869x -0.5755
Modelo exponencial es:
>> y1=exp(-0.5755)*x.^(1.8860)
>> X=linspace(0,5);
Y=-3.5100+2.9700*X;
x1=linspace(0,5);
>> y1=exp(-0.5755)*x1.^(1.8860);
>> plot(X,Y,'b',x1,y1,'g--'),grid
>> legend('Funcion Minimos Cuadrados','y=Cx^A')
Error
>> x=[1 2 3 4 5];
>> y=-3.5100+2.9700*x;
>> y1=exp(-0.5755)*x.^(1.8860);
>> Er=(1/5*((sum(y-y1)).^2)).^(1/2)
Er =
0.2262
El error en la curva ajustada de la función 𝑓(𝑥) = 𝐶𝑥 𝐴 presenta un error menor
con respecto al ajuste de mínimos cuadrados, se puede concluir que es la más
adecuada para estimar los valores
