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