Interpolación

2
A menudo se conocen los valores de la función f(x) en x0, x1, x2,...., xn, pero se
desconoce la expresión de f(x). La estimación de f(x), se conoce como
interpolación.

El objetivo del capítulo es estimar f(x) para cualquier punto x. Si x se encuentra

entre el valor menor y mayor de los xi se dice que el problema es de
interpolación; en cambio cuando x se encuentra fuera del rango, el problema
es de extrapolación.

Las clases de funciones aproximantes son: polinomios, exponenciales,

racionales, logaritmicas y trigonométricas.

En los años 1913 al 1923, Sir Edmund Wittaker desde la Universidad de

Edimburgo, desarrolló la teoría de interpolación en su cátedra de Matemática
Numérica.

MATLAB, tiene implementada las funciones que aparecen en la tabla

2.1 para realizar la interpolación.
.
Nombre Descripción

griddata Rejilla de datos

Interp1 Lookup tabla 1D
Interp2 Lookup tabla 2D
50

interpft Con método FFT

Tabla 2.1

El programa interpDemo.m de la figura 2.1, presenta el problema de

interpolación de una manera gráfica (figura 2.2); haciendo uso de la función
Interp1 mencionada en la tabla 2.1.

function InterpDemo()
% Demo Interpolacion
% E. Raffo Lecca
%
x=[0 1 4 6];
y=[1 -1 1 -1];
x0 = 0:0.10:6;
ylin = interp1(x,y,x0,'linear');
yspline = interp1(x,y,x0,'spline');

%ycubic = interp1(x,y,x0,'cubic');
plot(x,y,'o');
hold on;%para colocar varios plot;por default es off
plot(x0,ylin,x0,yspline);
legend('dato','lineal','spline',0);
%pos = 0,1,2,3,4,-1; % 0 = el mejor lugar o automático
title('Interpolación');

Figura 2.1: interDemo.m

2.1 Interpolación Polinómica

Los polinomios de interpolación, son los que más se usan en los cálculos por
computadora; y se basan en plantear un polinomio de grado n que pasa por los
puntos x0, x1, x2,..., xn. La razón se halla, en que las funciones son fáciles de
operar (las operaciones con el cálculo son simples) y proporcionan buenos
métodos, para la integración numérica, como a las ecuaciones diferenciales.

Los conjuntos de datos, como tablas de la forma:

 51

x0 x1 x2 … xn
y0 y1 y2 … yn

Se encuentran, envueltos en cálculos técnicos. La fuente de estos datos,

pueden ser observaciones experimentales, o de computación numérica. Se
asume implícitamente, que los puntos son precisos.

Existen dos métodos muy difundidos: el de Lagrange y de Newton.

Figura 2.2: Ejecución de interDemo.m

 52

Considere una familia de funciones de una variable x,

ϕ ( x ; a0 , a1 ,..., an ) ,

Donde existen (n+1) coeficientes: a0, a1, a2,..., an. El problema radica en
determinar los parámetros ai para los (n+1) pares de números (xi, fi), i = 1,...,
n con xi ≠ xk para i ≠ k , con la condición:
ϕ ( xi ; a0 , a1 ,..., an ) = f i , i = 0,..., n

Sea la forma lineal,

ϕ ( x ; a 0 , a1 ,..., a n ) = a 0ϕ0 ( x ) + a1ϕ1 ( x ) +... + a nϕn ( x )

Es aquí donde se encuentra la interpolación polinomial

ϕ ( x ; a0 , a1 ,..., an ) = a0 + a1 x + a2 x 2 + ... + an x n

2.1.1 Forma de Lagrange

Considere el par de puntos (x0, f(x0)) y (x1, f(x1)), aproximado mediante una
línea recta (ver figura 2.3); y constituye el caso mas simple de interpolación.

f(x1)

P(x)

f(x0)

x0 x x1
Figura 2.3: Interpolación lineal
 53

f ( x1 ) − f ( x0 ) P ( x) − f ( x0 )
=
x1 − x0 x − x0
( x − x0 )( f ( x1 ) − f ( x0 ))
P ( x) = f ( x0 ) +
( x1 − x0 )
x − x1 x − x0
P ( x) = f ( x0 ) + f ( x1 )
x0 − x1 x1 − x0

x1
x2
P(x)

x0
f(x) xn

x3

Figura 2.4: Polinomio de grado n=4

En la figura 2.4, se observa que para una interpolación de 5 puntos (n

vale 4); se obtiene un polinomio de grado 4.

Los (n+1) puntos de x, donde existe el valor de f (x), serán utilizados para
construir un polinomio P(x) de grado ≤ n que interpola a f (x) en los puntos x0,
x1, x2,..., xn, y satisface:

P(xi) = f (xi) , i = 0, 1, ... , n

Un polinomio de grado ≤ n que se anule en todos los puntos xi, excepto

en xk, resulta:
g k ( x ) = ∏ ( x − xi )
i =0
i≠k

= ( x − x0 )( x − x1 )...( x − xk −1 )( x − xk +1 )...( x − xn )
54

Además de:

g k ( xk ) = ∏ ( xk − xi )
i =0
i≠k

Luego se define un polinomio de Lagrange como:

g k ( x)
Lk ( x) =
g k ( xk )
x − xi
=∏
i =0 xk − xi
i≠k

Donde se observa que en la función cardinal Li ( x ) , se cumple:

1 i=j
Li ( x j ) = δij = 
0 i≠j

Donde, δij es el delta kronecker. Por ejemplo:

n n
P( x j ) = ∑ Li ( x j ) yi = ∑ yiδij = y j
i =0 i =0

El polinomio de grado ≤ n que interpola a f (x) en x0, x1, x2,..., xn es:

n
P ( x) = ∑ f ( xk ) Lk ( x)
k =0

P(x) se conoce como forma de Lagrange y Lk(x) son los polinomios de la

función cardinal, en los puntos x0, x1, x2,..., xn.

A modo de ejemplo usar la forma de Lagrange para producir un

polinomio cúbico y evaluar para x = 2 y 3.
 55

xi 0 1 4 6
yi 1 −1 1 −1

( x − x1 )( x − x2 )( x − x3 )
L0 ( x) =
( x0 − x1 )( x0 − x2 )( x0 − x3 )

( x − x0 )( x − x2 )( x − x3 )
L1 ( x) =
( x1 − x0 )( x1 − x2 )( x1 − x3 )

( x − x0 )( x − x1 )( x − x3 )
L2 ( x) =
( x2 − x0 )( x2 − x1 )( x2 − x3 )

( x − x0 )( x − x1 )( x − x2 )
L3 ( x ) =
( x3 − x0 )( x3 − x1 )( x3 − x2 )

Evaluando los polinomios de Lagrange para x = 2.

i f (xi) Li (x)

0 1 -0.33333
1 -1 1.06667
2 1 0.33333
3 -1 -0.06667

P(x) = -1

Para x = 3, se consigue P(3) = 0. Con la tabla:

i xi yi

0 1.05 1.02470
1 1.10 1.04881
2 1.15 1.07238
3 1.20 1.09544
56

El cálculo
de P(x), se presenta en la tabla:

I f (xi) Li (x)

0 1 -0.25000
1 -1 0.60000
2 1 0.75000
3 -1 -0.10000

Evaluando los polinomios de Lagrange para la tabla, se encuentra

P(1.12)= 1.0583.

i f (xi) Li (x)

0 1.02470 -0.064
1 1.04881 0.672
2 1.07238 0.448
3 1.09544 -0.056

Figura 2.5: Algoritmo de Lagrange

Figura 2.5: Algoritmo de Lagrange

 57

En la figura 2.5, se presenta el Algoritmo de Lagrange, y en la figura

2.6 el programa lagrange.m.

function p=Lagrange(x,y,x0)
% Interpolacion por Lagrange
% E. Raffo Lecca
% Datos
% x =es el vector x
% y =es el vector f(x)
% n =el orden de la matriz
% Resultados
% L =es el polinomio de lagrange
% p =valor interpolado
n = length(x);
L = zeros(n,1);
fprintf(' k x y L\n');
for k=1:n
prod1=1;prod2=1;
for i=1:n
if i~=k
prod1=prod1*(x0-x(i));
prod2=prod2*(x(k)-x(i));
end
end
L(k)=prod1/prod2;
end
p = 0;
for k=1:n
p = p+L(k)*y(k);
fprintf('%5d%10.6f%10.6f%10.6f\n',k,x(k),y(k),L(k));
end

Figura 2.6: Lagrange.m

El polinomio π( x ) de grado (n + 1), se define como el producto

π( x) = ( x − x0 )( x − x1 )( x − x2 )...( x − xn )

Ahora la derivada de π( x ) , es tal que uno de sus factores es

eliminado:
58

π' ( xi ) = ( xi − x0 )( xi − x1 )( xi − x2 )...( xi − xn )

Donde el factor en cuestión ( xi − xi ) , es omitido del producto.

La interpolación polinomial Lagrangeana, es de la forma

n
Pn ( x) = ∑Lk ( x) f ( xk )
k =0
n
π ( x)
=∑ f ( xk )
k =0 ( x − xk )π ' ( xk )

Donde:

π ( x)
Li ( x ) =
( x − xi )π ' ( xi )
( x − x0 )( x − x1 )...( x − xi −1 )( x − xi +1 )...( x − xn )
=
( xi − x0 )( xi − x1 )...( xi − xi −1 )( xi − xi +1 )...( xi − xn )

La interpolación polinomial es global; pero debe ser utilizada para un

número reducido de puntos. Para dos puntos cercanos, se utiliza la
interpolación lineal. Entre tres y seis puntos cercanos, se producen buenos
resultados. Más de seis puntos, se debe tomar con suspicacia. La razón, se
halla, en que la tendencia es oscilante.

Cuando las abscisas son equidistantes, los coeficientes Lagrangeanos,

han sido tabulados para varios valores del grado del polinomio(n).

Son muy usadas, aquellas donde n es un valor par (donde existen n +1

puntos). Para el caso de n + 1 = 2m + 1 , las abscisas son reenumeradas, como
x−m ,..., x−1 , x0 , x1 ,..., xm .

Siendo el espacio uniforme, con ancho h ; y éste es medido desde el

punto central; y siendo Li ( x) invariante a cualquier cambio lineal,
entonces:
 59

h = xi +1 − xi
x = x0 + hs , xi = x0 + hs i

( s − s0 )( s − s1 )...( s − si −1 )( s − si +1 )...( s − sn )
Li ( x) =
( si − s0 )( si − s1 )...( si − si −1 )( si − si +1 )...( si − sn )
( s + m)( s + m −1)...( s − i +1)( s − i −1)...( s − m +1)( s − m)
=
(i + m)( i + m −1)...( 2)(1)( −1)( −2)...( i − m +1)( i − m)

Entonces:

(1 − s 2 )(4 − s 2 )...(m 2 − s 2 )
L0 ( s) = ,y
(m!) 2

Li ( s) =
(− 1)i+1 s( s + i)
(m + i)!(m − i)!
[
(1 − s 2 )(4 − s 2 )...((i − 1) 2 − s 2 )((i + 1) 2 − s 2 )...(m 2 − s 2 ) ]

Para valores de i =±1, ±2,..., ±m . Siendo la nueva numeración

L−m ( x ), L−m+1 ( x ),..., Lm−1 ( x ), Lm ( x ) .

Los coeficientes de Lagrange para un polinomio de grado 3, es

conocido como Lagrange de 3 puntos; y para valores de m = 1, se presenta en
la tabla 3.2.

Esta tabla considera, el hecho que:

Li ( −s ) = L−i ( s )

s L−1 ( s ) L0 ( s ) L1 ( s )
0.0 0 1 0 0.0
0.1 -0.045 0.99 0.055 -0.1
0.2 -0.080 0.96 0.120 -0.2
0.3 -0.105 0.91 0.195 -0.3
0.4 -0.120 0.84 0.280 -0.4
0.5 -0.125 0.75 0.375 -0.5
 60

0.6 -0.120 0.64 0.480 -0.6

0.7 -0.105 0.51 0.595 -0.7
0.8 -0.080 0.36 0.720 -0.8
0.9 -0.045 0.19 0.855 -0.9
1.0 0 0 1 -1.0
L1 ( s ) L0 ( s ) L−1 ( s ) s

Tabla 3.2: Coeficientes de Lagrange para 3 puntos

Desde la tabla, interpolar con Lagrange a tres decimales, cuando x =

1.164 y 1.142.

x 1.14 1.15 1.16 1.17

y 0.1185 0.1189 0.1213 0.1230

a) x = 1.164, es cercano a 1.16; y

h =1.17 −1.16
De x = a + sh , con a =1.16
1.164 −1.16
s= = 0.4
0.01

Pn(x) = -0.12 (0.1189) + 0.84 (0.1213) + 0.28 (0.1230) = 0.1221

b) x = 1.142, es cercano a 1.14; y

h =1.16 −1.15
De x = a + sh , con a =1.15
1.142 −1.15
s= = −0.8
0.01

Pn(x) = 0.72 (0.1185) + 0.36 (0.1189) -0.08 (0.1213) = 0.1184

2.1.2 Forma de Newton

 61

Sea un polinomio de grado n-ésimo:

Pn ( x) = b0 + b1 ( x − x0 ) + b2 ( x − x0 )( x − x1 ) + ...
... + bn ( x − x0 )( x − x1 )....( x − xn −1 )

Haciendo Pn (xi) = fi , i = 0, 1,..., n

f 0 = Pn ( x0 ) = b0
f1 = Pn ( x1 ) = b0 + b1 ( x − x0 )
f 2 = Pn ( x2 ) = b0 + b1 ( x − x0 ) + b2 ( x − x0 )( x − x1 )
...

Luego se definen:
f1 − f 0
f 01 =
x1 − x0
f 2 − f1
f12 =
x2 − x1
f12 − f 01
f 012 =
x2 − x0

El polinomio es equivalente a:

Pn ( x) = f 0 + f 01 ( x − x0 ) + f 012 ( x − x0 )( x − x1 ) + ...
... + f 012 ... k ( x − x0 )( x − x1 )...( x − xk −1 )

De f 012...k = f [x0, x1,..., xk], se tiene que los coeficientes son:

b0 = f ( x0 )
b1 = f [ x0 , x1 ]
b2 = f [ x0 , x1 , x2 ]
...
bn = f [ x0 , x1 ,..., xn ]
62

Donde f01...k se conoce como diferencia dividida. La primera diferencia

dividida se calcula como:

f ( x j ) − f ( xi )
f [ xi , x j ] =
x j − xi

La segunda diferencia dividida, es la diferencia de las dos primeras diferencias

divididas

f [ x j , xk ] − f [ xi , x j ]
f [ xi , x j , xk ] =
xk − xi

Y en general:
f [ x1 , x2 ,..., xk ] − f [ x0 , x1 ,..., xk −1 ]
f [ x0 , x1 , x2 ,..., xk ] =
xk − x0

En la tabla 2.2, se presenta el esquema para las diferencias divididas

es; que permite una gran facilidad de cálculo.

k =1 k =2 k
x0 f0 f 01 f 012 f 0123
x1 f1 f12 f123
x2 f2 f 23
... ...
xn fn

Tabla 2.2: Diferencias divididas

En la figura 2.7, se presenta el programa InterpolNewton.m, que

ejecuta la interpolación de Newton.

Sean los siguientes valores:

xi 1 4 6 5
yi 0 1.3862944 1.7917595 1.6094379
 63

f ( x1 ) − f ( x0 ) 1.3862944
f [ x0 , x1 ] = = = 0.46209813
x1 − x0 4 −1

1.7917595 − 1.3862944
f [ x1 , x2 ] = = 0.20273255
6−4

1.6094379 − 1.7917595
f [ x2 , x3 ] = = 0.18232160
5−6

f [ x1 , x2 ] − f [ x0 , x1 ]
f [ x0 , x1 , x2 ] = = −0.05187312
6 −1

f [ x2 , x3 ] − f [ x1 , x2 ]
f [ x1 , x2 , x3 ] = = −0.02041095
5−4

f [ x1 , x2 , x3 ] − f [ x0 , x1 , x2 ]
f [ x0 , x1 , x2 , x3 ] = = 0.00786554 ,
5 −1

La tabla de diferencias divididas finitas es:

i xi f (xi) 1 2 3

0 1 0 0.46209813 -0.05187312 0.00786554

1 4 1.3862944 0.20273255 -0.02041095

2 6 1.7917595 0.18232160
3 5 1.6094379

La forma de Newton es:

Pn ( x ) = 0 + 0.46209813 ( x −1) − 0.05187312 ( x −1)( x − 4) +

0.00786554 ( x −1)( x − 4)( x − 6)
Pn ( 2) = 0.62877

Sean los valores:

64

i 0 1 2 3
xi 0 1 4 6
yi 1 −1 1 −1

La tabla de diferencias divididas finitas es:

i xi f (xi) 1 2 3

0 0 1 -2.0 0.66666667 -0.16666667

1 1 -1 0.66666667 -0.33333333
2 4 1 -1.0
3 6 -1

P ( x ) =1.0 −2( x −0) +0.66666667 ( x −1) x −....

−0.16666667 ( x −1)( x −4) x
P ( 2) =−1.0

function p=InterpolNewton(x,y,x0)
% Interpolacion por Newton
% E. Raffo Lecca
%
n=length(x);
T=zeros(n,n);
for i=1:n
T(i,1)=y(i);
end
%T(:,1)=y';
for j=2:n
for i=1:n+1-j
T(i,j)=(T(i+1,j-1)-T(i,j-1))/(x(i+j-1)-x(i));
end
end
%Evaluar P(x)
suma=0;
producto=1;
for j=1:n
suma=suma+T(1,j)*producto;
producto=producto*(x0-x(j));
 65

end
p=suma;
%impresion
fprintf(' X Difer. Divididas\n');
fprintf('==========');
for j=1:n
fprintf('==========');
end
fprintf('\n');
for i=1:n
fprintf('%10.3f',x(i));
for j=1:n+1-i
fprintf('%10.3f',T(i,j));
end
fprintf('\n');
end

Figura 2.7: InterpolNewton.m

Al aplicar, la regla de Horner, para evaluar el polinomio de

interpolación de Newton, se encuentra:

cn = bn
ck = bk + ( x-x k )bk +1 , k = n −1, n − 2,..., 0

Para P ( 2) = −1.0 , los cálculos, son como aparece a continuación:

k xk bk ck
3 6 -0.16666667 -0.16666667
2 4 0.66666667 1
1 1 -2.0 -1
0 0 1 -1

En el problema propuesto 15 del capítulo, se presenta un ejemplo de

interpolación inversa.

2.1.3 Algoritmo de Neville

66

Una forma de resolver el problema de interpolación, es considerar resolver el

problema para un pequeño conjunto de valores; y después actualizar estas
soluciones para obtener la solución completa.

Sea por el conjunto de puntos (xi , fi), i = 0,1,..., n; se denota por

Pi0 i1 ... ik al polinomio, donde

Pi0 i1 ...ik ( xi j ) = f i j , j = 0,1,...., k

Por recurrencia, se tiene:

Pi ( x ) = f i
( x − x i0 ) Pi1 ... ik ( x ) − ( x − x ik ) Pi0 i1 ... ik −1 ( x )
Pi0 i1 ... ik ( x ) =
xik − x i0

k =0 k =1 2 3
x0 f 0 = P0 ( x)
P01 ( x )
x1 f1 = P1 ( x) P012 ( x )
P12 ( x) P0123 ( x)
x2 f 2 = P2 ( x) P123 ( x)
P23 ( x )
x3 f 3 = P3 ( x)

El cálculo para P012(x), viene dado por

( x − x0 ) P12 ( x) − ( x − x2 ) P01 ( x)
P012 ( x) =
x2 − x0

Determinar P0123(2):
 67

k =0 1 2 3
0 1
−3
1 −1 1
1 1
2 1 1
1
3 −1

( 2 − 0)(−1) − (2 − 1)(1)
P01 (2) = = −3
1− 0

(2 − 1)(1) − (2 − 2)(−1)
P12 (2) = =1
2 −1

(2 − 2)(−1) − (2 − 3)(1)
P23 (2) = =1
3− 2

(2 − 0)(1) − ( 2 − 2)(−3)
P012 (2) = =1
2−0

(2 − 1)(1) − (2 − 3)(1)
P123 (2) = =1
3 −1

(2 − 0)(1) − (2 − 3)(1)
P0123 (2) = =1
3−0

Una manera de llevar el cálculo para Neville, de una manera cómoda,

es seguir la siguiente receta:

Ti + k , k = Pi , i +1,..., i + k
68

Con tiene la tabla

x0 f 0 = T00
T11
x1 f1 = T10 T22
T21 T33
x2 f 2 = T20 T32
T31
x3 f 3 = T30

Donde:

Ti 0 = f i

( x − xi − k ) Ti , k −1 − ( x − xi ) Ti −1, k −1
Tik =
( xi − xi − k )

Ti , k −1 − Ti −1, k −1
= Ti , k −1 + , 1 ≤ k ≤ i , i = 0,1,..., n
x − xi −k
−1
x − xi

2.1.4 Error en una interpolación polinomial

Sea la función f (x); con sus valores

f i = f ( xi ) , i = 0, 1, ..., n

Los cuales son interpolados. Una interpolación polinomial P(x) =P0...n(x) con

P ( xi ) = f i , i = 0, 1, ..., n

Reproduce f (x), para el argumento xi.

 69

El error es medido por f ( x ) − P( x) ; p

ara valores x ≠ xi ,i= 0,1,..., n.

Teorema:
Si una función f tiene (n+1) derivadas; entonces para cada argumento X ,
existe un número
ξ en el más pequeño intervalo I [xo, ...., xn, X ] que contiene a X y a todas
las abscisas xi, satisfaciendo:

π( x) f ( n +1)
(ξ )
f ( x ) − P01 ... n ( x ) =
( n +1)!

Donde π(x ) también conocido como polinomio soporte, es un

polinomio de grado n. Soporte, se conoce a las abscisas xi.

π( x ) = ( x − x0 )( x − x1 )....( x − xn )

Así:
f ( x) =cos x
iπ
xi = , i = 0,1,2,3,4,5 o n =5
10
cos ξ
f ( x) − P ( x ) = ( x − x0 )( x − x1 )...( x − x5 )
720
| π( x ) |
f ( x) −P ( x ) ≤ , donde ξ ∈I [ x0 ,...., xn , x)
720

Por la fórmula fundamental de Newton, el valor real es igual al valor

interpolado más el resto o residuo:

f ( x ) = Pn ( x) + Rn ( x )

Siendo Pn (x ) el polinomio de interpolación de grado n.

Pn ( x) = f [ x0 ] + ( x − x 0 ) f [ x 0 x1 ] + ( x − x 0 )( x − x1 ) f [ x0 x1 x 2 ] + ...
+ ( x − x 0 )( x − x1 )...( x − x n −1 ) f [ x 0 x1 x 2 ... x n ]
70

Y residuo
Rn ( x ) = ( x − x0 )( x − x1 )...( x − x n −1 )( x − xn ) f [ x0 x1 x2 ... xn x n +1 ] ;
para el caso de existir el punto ( x n +1 , f ( x n +1 )) .

n
Rn ( x) = ∏[( x − xi )] f [ x0 x1 x 2 ... x n x n +1 ]
i =0

Teorema de Rolle:

Si f(x) es continua para a ≤ x ≤ b y diferenciable para a < x < b ; si f(a)

=f(b) entonces f’( ξ ) = 0 para al menos un ξ en a < ξ < b.

Luego:
f ( x ) = Pn ( x ) + Rn ( x ), y
Q (t ) = f (t ) − Pn (t ) − Rn ( x )

Se observa que para t = xi , i = 0,1,2,..., n, Q ' (t ) = 0 ; por ser

f ( xi ) = Pn ( xi ); siendo Rn ( x) = 0 para esos puntos.

Definiendo el resto para cualquier punto como

Rn ( x) = π( x ) K
Entonces:

 n 
Q (t ) = f (t ) − Pn (t ) − ∏ ( x − xi ) K
 i =0 

Al aplicar el teorema de Rolle a Q(t), se obtiene Q ' (t ) = 0 ; y como

por definición Pn(t) es un polinomio de grado n, entonces su derivada (n + 1)
es igual a cero:

Q ( n+1) (t ) = f ( n +1)
(t ) − P ( n+1) (t ) −( n +1)! K

de donde:
 71

( n +1)
f (t ) −( n +1)! K = 0
f ( n+1) (t )
K =
( n +1)!

y el residuo aparece en una fórmula más compacta:

Rn ( x ) = π ( x ) K
f ( n+1) (t )
= π ( x)
(n +1)!

Para algún t = ξ en el intervalo.

Teorema:
Si xo,..., xn son (n+1) puntos distintos contenidos en el intervalo
[a,b] y f (ξ) es una función derivable n+1 veces en (a, b), entonces:

π( x ) f ( n +1)
(ξ )
f ( x) − P ( x) ∞ ≤
( n +1)!

El polinomio soporte, para formalizar su notación, se expresará como

πn (x ) , para asociarlo con Pn ( x ) que es de grado n.

Para puntos equidistantes, con x = a + th , se expresa como:

πn ( x) = h n+1t (t −1)( t − 2)...( t − n)

Siendo su error máximo:

π n ( x) ∞[ 0,n ] = h n+1 t (t −1)(t − 2)...( t − n) ∞[ 0 ,n ]

Por ejemplo:
72

h2
π1 ( x) ∞[ 0,1] = h 2 t (t −1) ∞ =
[ 0 ,1 ] 4
2h 3
π2 ( x) ∞[ 0, 2 ] = h 3 t (t −1)( t − 2) ∞[ 0, 2 ] =
3 3
π3 ( x) ∞[ 0,3] = h t (t −1)( t − 2)( t − 3) ∞[ 0,3] = h 4
4

El polinomio soporte, es posible estimarlo, en una cota superior de

n!
πn ( x) ∞[ a ,b ] ≤ h ( n+1) ; teniendo los casos particulares de n =1, 2, 3.
4

Cuando el error, se está midiendo en un punto, entre el mínimo y el

máximo en la subdivisión, se puede aplicar la desigualdad:

( n +1)
f
f ( x ) −P ( x ) ≤ ∞
h ( n+1)
∞
4( n +1)
1
Por ejemplo, para f ( x) = , aproximado a una cuadrática:
2 x
1
f ( x) = , con x ∈[θ ,1], y n = 2
2 x
15 1 −θ
f (ξ) = − ξ −7 / 2 , h =
( 3)

16 2
(1 − θ ) 3
15
f ( x) − P ( x ) ∞ ≤
2 .3 16 θ 7 / 2
5

2.2 Polinomios de Chebyshev

La ortogonalidad, es una generalización del concepto, de la perpendicularidad,
a un espacio vectorial, con producto interno.

El producto interno, en un espacio vectorial V, es una operación en V,

que asigna a un par de vectores x e y en V, un número real (x, y), satisfaciendo
ciertas condiciones.

En el caso del espacio vectorial Rn, es el producto escalar (x, y)=xT y.

 73

2.2.1 Polinomios ortogonales

Sea un conjunto de polinomios p1 ( x ), p2 ( x ),..., pn ( x ) , donde cada
miembro es ortogonal a todos los otros miembros en el conjunto.

Este conjunto, tiene como miembros a

∏ ={ pj n ( x ), n ≤ j }

y su producto interno en C[a, b] (el espacio vectorial de funciones definidas y

continuas), se define como:

b
( f , g ) = ∫ w( x ) f ( x ) g ( x ) dx
a

Tales polinomios cumplen con la propiedad ortogonal

( pi , p j ) = 0, si i ≠ j

Siendo únicamente definidos por recurrencia.

Si el polinomio pn(x) sobre (a, b), es ortogonal a todos los polinomios

de grado inferior a n; entonces se cumple que el producto escalar (p, q) es
cero.

b
( p, q ) = ∫ w( x ) pn ( x) qn −1 ( x) dx = 0
a

Donde el sistema de polinomios ortogonales, está asociado con

w( x ) , una función de peso; con qn −1 ( x ) , un polinomio de grado menor o
igual que n -1.

El término aproximación de Chebyshev, es muy utilizado para

expresar la minimización del error máximo; es decir el criterio del minimax.
74

2.2.2 Aproximación minimax

En la búsqueda de los nodos x0 , x1 , x2 ,..., xn ∈[ a, b ] , para que la
interpolación polinomial Pn sea el mejor polinomio de aproximación para una
función f(x) sobre el intervalo [a, b]; tal aproximación es:

min max f ( x ) − Pn ( x)
p a ≤x ≤b

Conocida como la aproximación minimax.

El error de la interpolación por polinomios, viene dado por:

π( x) f ( n +1)
(ξ )
Rn ( x ) =
( n +1)!

Donde ξ es el punto de valor medio en el intervalo [a, b] y

n
π ( x ) = ∏ ( x − xk ) .
k =0

En la consideración de ser f un polinomio mónico (un polinomio

donde an , el coeficiente líder, vale 1), de grado (n+1), es decir:

f ( x) = a0 + a1 x + a2 x 2 + ... + x n +1 ,

Se tiene que f (n+1)(x) = (n+1)! y el error se reduce a Rn ( x ) = π( x) . Esto

significa que la búsqueda de los nodos x0 , x1 , x2 ,..., xn se simplifica,
cuando max π( x) se hace mínimo.

La solución a este problema, queda a cargo de los polinomios de

Chebyshev o polinomios Tn. Pafnuti Chebyshev (1821-1894), fue un
matemático ruso.

2.2.3 Polinomios Tn
En Chebyshev para el intervalo (-1,1), la función de peso es:
 75

1
w( x) =
1 − x2

El polinomio pn (x ) , es tal que,

1
pn ( x ) qn −1 ( x )
∫
−1
1 − x2
dx = 0

Haciendo el cambio de variables x = cos θ , se tiene

∫p
0
n (cos θ ) qn −1 (cos θ ) dθ = 0

Si qn −1 (cos θ) se puede expresar como una combinación

1, cos θ,..., cos nθ , entonces

∫p
0
n (cos θ )(cos kθ ) dθ = 0, k = 0,1,..., n −1

De donde se deduce que:

pn ( x) = A cos nθ

Cuando A = 1, estos polinomios son conocidos como de Chebyshev(o

Tschebyscheff); y se escriben como:

pn ( x ) = cos( n cos −1 x) = Tn ( x )

El producto escalar tiene la propiedad ortogonal:

1
Tr ( x )Tk ( x )
∫
−1 1 − x2
dx = 0, r ≠ k
76

Un polinomio es una combinación lineal de monomios

1, x, x 2 ,..., x n . Cada monomio en el intervalo [-1,1], tiene un máximo con
magnitud de 1 en x = ±1 y un mínimo igual a cero en x = 0.

Sea el polinomio Pn(x), que aproxima a una función f(x), y definido

por:

Pn ( x) = a0 + a1 x + a2 x 2 + ... + an x n

Donde los coeficientes, son a0 , a1 , a2 ,..., an ; y x se encuentra en [-1,1].

Para la situación de f(z), con z sobre el intervalo [a,b], la

transformación al intervalo − 1 ≤ x ≤ 1 , es por el cambio lineal de variable:

b −a z −a
=
1 − (−1) x − ( −1)
2z −b − a
x= , z ∈[ a, b] ,
b −a

Esta transformación permite examinar los monomios sobre el intervalo

[-1,1]. La aproximación producida por el polinomio permitirá reducir el error
máximo a un valor mínimo.
 77

Figura 2.8: Coseno de nθ

Examinando la figura 2.8, para las funciones cosθ, cos2θ, se observa

que tienen los mismos valores de mínimos y máximos en el intervalo
0 ≤ θ ≤ π . Por una simple transformación cos nθ en [0,π] es Pn(x) en [-1,1].

Los polinomios de Chebyshev de primera clase, son definidos como:

T0 =1
Tn ( x ) = cos( nθ ), n =1,2,...

De la figura 2.9, se desprende que x = cos θ y θ = arccos(x). Luego:

78

cos θ = cos(arccos x) = x
T1 ( x) = x
cos 2θ = 2 cos 2 θ −1
T2 ( x) = 2( cos θ ) −1
2

= 2 x 2 −1

Y en general la identidad trigonométrica:

cos nθ = 2 cos θ cos( n −1)θ − cos( n − 2)θ ,

Es usada para calcular los polinomios de Chebyshev de alto orden; dando

origen a la siguiente relación de recurrencia:
Tn ( x ) = 2 xT n −1 ( x ) −Tn −2 ( x )

-1 1 x
Figura 2.9

El cálculo para T3(x), es como sigue:

T3 ( x ) = 2 xT 2 ( x ) −T1 ( x )
= 2 x ( 2 x 2 −1) − x
= 4 x 3 −3x

Tn Expresión
T0 1
T1 x
T2 2x2 - 1
T3 4x3 - 3x
T4 8x4 - 8x2 + 1
T5 16x5 - 20x3 +5x
 79

T6 32x6 – 48x4+ 18x2 – 1

Tabla 2.3: Polinomios de Chebyshev

Figura 2.10

En la tabla 2.3, se presentan otros polinomios de Chebyshev; y en la

figura 2.10, las formas geométricas mas comunes.

Sea Tn un polinomio de grado n; su coeficiente an, o líder, tiene como

valor a 2n-1, n≥1.

Teorema:

max 2 −n Tn +1 ( x ) = 2 −n . max Pn +1 ( x )
−1≤x ≤1 −1≤x ≤1

Pn+1(x) es un polinomio grado n+1. El polinomio Tn+1(x)2-n es un polinomio

mónico. El polinomio Tn+1(x) tiene un máximo con valor absoluto de 1; porque
cos[( n +1) arccos x ] = ±1 : Luego se cumple:
80

2 −n Tn+1 ( x ) =2 −n Tn+1 ( x ) ∞
=2 −n
∞

Desde el teorema anterior, se desprende que Tn+1(x) es el polinomio

que se está buscando, y posee (n+1) raíces, los soportes S =
{ x0 , x1 , x2 ,..., xn } . Estos nodos son los que utilizan las interpolaciones de
Lagrange o la diferencia dividida de Newton para formar un polinomio de
grado n.

Se cumple que:

max π( x) ≤ max Pn +1 ( x)
−1≤x ≤1 −1≤x ≤1

La relación T 'n (1) = n , se cumple en el intervalo [-1,1].

2

Figura 2.11

El lema de propiedad minimax, plantea, que de todos los polinomios

1−n
mónicos, el polinomio 2 Tn ( x) , tiene la magnitud más pequeña en [-
1,1]. En la figura 2.11, se presenta el polinomio soporte, para puntos
equidistantes, y para los puntos de Chebyshev, en donde se observa, grandes
 81

oscilaciones, para el polinomio soporte, con puntos equidistantes. Se ha usado

n = 10. Se concluye que para los puntos Tn ( x ) ,

πn ( x) ∞[ −1,1] = 2 −n Tn+1 ( x) ∞ = 2 −n .

Teorema minimización de soporte

Sean Tn+1 ( x ) el polinomio de Chebyshev de grado n+1 y el soporte
formado por sus raíces S = { x0 , x1 , x2 ,..., xn } ; entonces se tiene que:

πS ( x ) ∞[ −1,1] ≤ πS ( x) ∞[ −1,1] , ∀S = { x0 , x1 ,... xn } ⊂[−1,1]

Donde S el soporte formado por los puntos equidistantes.

Figura 2.12

El conjunto de n puntos o raíces de los polinomios de Chebyshev, son

calculados de acuerdo a la siguiente relación:

( 2k −1)π
xk = cos , k =1,..., n
2n

Esta relación obedece a la siguiente deducción basada en un polinomio

de grado n:

Tn ( x ) = cos( nθ) = cos[ n(cos −1

( x))], de x = cos θ
Tn ( x ) = cos[ n cos −1 ( x )] = 0, para sus raices
n cos −1 ( x ) = cos −1 (0) ={π / 2,3π / 2,5π / 2,...}
n cos −1 ( x ) = ( 2k +1)π / 2, k = 0,1,..., n −1
( 2k +1)π
xk = cos , k = 0,1,..., n −1
2n

Por ejemplo para T3(x), las raíces desde 4 x 3 − 3 x = x (4 x 2 − 3) son:

x = 0 y x = ± 3 / 2 = ± 0.8660. Ver figura 2.12.
82

Equivalentes a: cos π/6, cos π/2 y cos 5π/6. Para la situación de

plantear un polinomio de grado a lo sumo 2 que sea un interpolante de
f ( x ) = x 3 + x +1 ; se tiene que teniendo los pares (x,y) y la aplicación de las
diferencias de Newton, se consigue el polinomio de grado 1. Como se presenta
en la figura 2.13.

Figura 2.13: Usando interpolación de Newton

El polinomio interpolante es: P2 = -0.516 + 1.750(x-x0) = 1 + 7x/4.

Usar los ceros de T3, para construir un polinomio interpolante de grado

dos, para f ( x) = ln x, en [1,2].

En la figura 2.14, se presenta la solución usando diferencias divididas,

a los nodos de T3(x). Como el intervalo es [1,2], o en términos generalizados es
[a, b], se hace necesario de la transformación lineal:

− 1
xk = [ (b − a) x k + a + b]
2
 83

Figura 2.14: aplicando la interpolación de Chebyshev

El polinomio resultante, P2(x) = 0.659080 + 0.585698(x-x0)-

0.232030(x-x0)(x-x1). En la figura 2.15, se presenta el M-file Chebyshev.m.

function p=InterpolChebyshev(f,n,a,b,x0)
% Interpolacion por Chebyshev
% E. Raffo Lecca
%
n=n+1;
x=zeros(n,1);
y=zeros(n,1);
for k=1:n
i=k-1;
x(k)=cos(0.5*(2*i+1)*pi/n);
end
x=0.5*((b-a)*x+a+b);
for k=1:n
y(k)=f(x(k));
end
InterpolNewton(x,y,x0);

Figura 2.15: M-file Chebyshev.m

(En este momento, refiérase a los problemas del capítulo).

2.2.4 Ortogonalidad
Dos miembros de una función, se dice que son ortogonales con respecto a una
función de peso w, si su producto interno es cero.

Examinando las propiedades de ortogonalidad:

π 0 ( m ≠ n)

∫0 cos mθ cos nθ dθ = π / 2 (m = n ≠ 0)
π ( m = n = 0)


Para el caso de los polinomios de Chebyshev, se tiene:

84

1 0 ( m ≠ n)
Tm ( x )Tn ( x) 
∫ 1−x2
dx = π / 2 ( m = n ≠ 0)
−1 π ( m = n = 0)


Esto es, los polinomios de Chebyshev, forman un conjunto ortogonal sobre el

intervalo -1≤ x ≤ 1, con una función de peso definida por:

1
w( x ) =
1− x2

Figura 2.16: La función de peso

La función de peso w(x) para los polinomios de Chebyshev, tiene la

particularidad, que incrementa la ponderación relativa, cuando los puntos se
encuentran en los límites del intervalo. En la figura 2.16, para el intervalo [-
1,1], se observa que: w(0) = 1 y w(x) → ∞, cuando x = ± 1.
 85

2.3 Interpolación de Hermite

La interpolación polinomial es global; es decir se usa una función polinomial
para pasar a través de todos los datos. Cuando se adicionan nuevos puntos, se
requiere de incrementar el grado del polinomio.

Desde 1960, un método alternativo se ha hecho muy popular: las

funciones polinómicas por pieza (piecewise polynomial functions). Los spline
cúbicos y los hermite son ejemplos de estas funciones.
Los spline son usados para resolver ecuaciones diferenciales.

2.3.1 Interpolación lineal de Hermite

Una función polinomial lineal por pieza definida para las x, tiene la propiedad
que L(x) es una línea quebrada entre xi y xi+1 ; vale decir que L(x) permite
diferentes líneas entre cada par de juntas adyacentes (del inglés Knots, joints o
breakpoints).

La propiedad
L( xi ) = yi , i = 0,1,2,..., n

Indica que

x − xi +1 x − xi , xi ≤ x ≤ xi +1
L ( x ) = yi + yi +1
xi − xi +1 xi +1 − xi , i = 0,1,2,...., n − 1

Un polinomio de interpolación por pieza, tiene por propiedad que si

los valores de yi son conocidos desde una función conocida g(x), se pueden
adicionar más puntos entre x0 y xn; y el interpolante converge a la función
original.

En la figura 2.17, se presenta la función HermiteLin.m para el Hermite

lineal. En la figura 2.18, se presenta a HermiteLinDemo.m, para efectuar un
demo de HermiteLin.m.
86

Si los datos yi son valores de la función g(x), y tiene una segunda

derivada continua, se demuestra que:

| L( x) − g ( x) | < 18 h 2 m á x| g´´(x) |= O(h 2 ) .

function L = HermiteLin(x,y,x0)
% Hermite lineal
% E. Raffo Lecca
% Datos
% x =es el vector x
% y =es el vector f(x)
% n =el orden de la matriz
% Resultados
% L =valor interpolado
n=length(x);
if x0<=x(1)
L=y(1);
else
if x0<=x(n)
i=1;
while x0>x(i)
i=i+1;
end
j=i-1;
L=y(j)*(x0-x(i))/(x(j)-x(i))+y(i)*(x0-,... x(j))/
(x(i)-x(j));
else
L=y(n);
end
end

Figura 2.17: HermiteLin.m

En la figura 2.19, se presenta el polinomio interpolante para los puntos

(x,y):
x 0 1 4 6
y 1 −1 1 −1
 87

function HermiteLinDemo()
% Demo de Interpolacion Hermite Lineal
% E. Raffo Lecca
% Datos
% x = es el vector x
% y = es el vector f(x)
% xinterpol = valor de x a interpolar
% Resultados
% yinterpol = valor interpolado
x=[0 1 4 6];y=[1 -1 1 -1];
xinterpol=0:0.1:6;
%invocacion a Hermite Lineal
for i=1:length(xinterpol)
yinterpol(i)=HermiteLin(x,y,xinterpol(i));
end
%grafica
plot(x,y,'o',xinterpol,yinterpol);
grid on;
xlabel('X');ylabel('Y');
title('comparacion con la interpol de Hermite Lineal');
legend('real','interpol Hermite Lineal',2);

Figura 2.18: HermiteLinDemo.m

88

Figura 2.19: Ejecución de HermiteLinDemo.m

2.3.2 Interpolación cúbica de Hermite

Sea una función C(x), formada por n segmentos de polinomios de grado 3:
p0(x), p1(x), p2(x),..., pn(x); con

pi ( x ) y x ∈[ xi , xi +1 ] ∀i = 0,1,..., n −1

Definida por:
 89

 p0( x) y x ∈ [ x0, x1]

 p1( x) y x ∈ [ x1, x 2]

C ( x) = 
...
 pn − 1( x) y x ∈ [ xn − 1, xn]

Por concepto de interpolación polinómica

C ( xi ) = f ( xi ), i = 0,1,2,..., n
De donde se deduce que:

pi ( xi ) = f ( xi ), i = 0,1, 2,..., n −1
pi ( xi + 1) = f ( xi + 1), i = 0,1,2,..., n −1

1 2 3
xi fi f[xi xi] f[xixi xi+1] f[xixixi+1xi+1]
xi fi f[xi xi+1] f[xixi+1 xi+1]
xi+1 fi+1 f[xi+1 xi+1]
xi+1 fi+1

Para evaluar f[xixi], se hace uso del concepto de la derivada como una
aproximación; así para ∆xi muy pequeño, se tiene:

f ( xi + ∆xi ) − f ( xi )
f [ xi x i] = lim Δx →0 = f ' ( xi )
∆xi

Luego

b 0 = f ( xi )
b1 = f ' ( xi )
f [ xixi + 1] − b1
b2 =
∆xi
f ' ( xi + 1) − 2 f [ xixi + 1] + b1
b3 =
∆x 2 i
90

En las figuras 2.20 y 2.21 se presentan los programas MATLAB para

la interpolación de Hermite cúbica; y en 2.22 una ejecución del mismo.

function valor=Hermite(x,y,yderiv,x0)
% Interpolacion segementaria por Hermite
% E. Raffo Lecca
%
% Datos
% x =es el vector x
% y =es el vector f(x)
% yderiv=es el vector de derivadas
% x0 =valor de x a interpolar
% n =el numero de segmentos
% np =numero de polinomios cubicos p
% Resultados
% C =es el conjunto de todos los polinomios p
% valor =valor interpolado
np=length(x);
n=np-1;
for i=1:n
C(1,i)=y(i);
C(2,i)=yderiv(i);
dx=x(i+1)-x(i);
difer=(y(i+1)-y(i))/dx;
C(4,i)=(yderiv(i+1)-2*difer+C(2,i))/dx/dx;
C(3,i)=(difer-C(2,i))/dx-C(4,i)*dx;
end
Ok=0;
if x0<x(1)
Ok=1;
j=1;
else
for i=2:np
if x0<x(i)
Ok=1;
j=i-1;
break;
end
end
end
if Ok==0
j=n;
end
difer = x0-x(j);
 91

valor = C(1,j)+C(2,j)*difer+C(3,j)*difer^2+C(4,j)*(x0-
x(j+1))*difer^2;

Figura 2.20: Hermite.m

function HermiteDemo()
% Demo de Interpolacion segementaria por Hermite
% E. Raffo Lecca
% Datos
% x =es el vector x
% y =es el vector f(x)
% xderiv =es el vector de derivadas
% xinterpol =valor de x a interpolar
% Resultados
% yinterpol =valor interpolado
x=0:pi/8:pi/2;y=cos(x);
xderiv=-sin(x);
xinterpol=0:pi/50:pi/2;
%invocacion a Hermite
for i=1:length(xinterpol)
yinterpol(i)=Hermite(x,y,xderiv,xinterpol(i));
end
%grafica
plot(x,y,'o',xinterpol,yinterpol);
grid on;
xlabel('angulo rad');
ylabel('cos');
title('comparacion cos y la interpolacion de Hermite');
legend('cos','interpol Hermite')

Figura 2.21: HermiteDemo.m

92

Figura 2.22: Ejecución de HermiteDemo.m

2.4 Interpolación segmentaria o spline

Cuando se interpola polinomios de alto orden, el tamaño del error también es
considerable. Se sugiere dividir el intervalo de los puntos x0, x1,..., xn en
varios segmentos y buscar el polinomio interpolante para cada caso; haciendo
que los cambios entre ellos no sean abruptos, sino suavizados.

La interpolación segmentaria o spline tiene su origen en las plantillas que

hacen uso los dibujantes para suavizar las curvas entre los puntos de un
segmento.
 93

2.4.1 Interpolación segmentaria lineal

Definiendo:
f ( xi +1 ) − f ( xi )
mi = ,y
xi +1 − xi

f ( x ) − f ( xi )
mi =
x − xi
Se tiene la forma
f ( x ) = f ( xi ) + mi ( x − xi )

f ( x ) = f ( x 0 ) + m0 ( x − x 0 ) x 0 ≤ x ≤ x1
f ( x ) = f ( x1 ) + m1 ( x − x1 ) x1 ≤ x ≤ x 2

f ( x ) = f ( x n −1 ) + mn −1 ( x − x n −1 ) x n −1 ≤ x ≤ x n

x 3.0 4.5 7.0 9.0

f (x ) 2.5 1.0 2.5 0.5

En x = 8.0 (ver figura 2.23)

f ( x ) = f ( x2 ) + m2 ( x − x2 ) x2 ≤ x ≤ x3

f ( x3 ) − f ( x 2 ) 0.5 − 2.5
m2 = = = −1
x3 − x 2 9.0 − 7.0

f (8) = f ( x 2 ) −1(8 − x 2 ) = 2.5 −1(8 − 7) = 1.5

La expresión general para un spline lineal es:

S(x) = yi + mi (x – xi) para xi ≤ x ≤ xi+1

94

Donde mi = ( y i+1 – y i) / (x i+1 – x i) , y

 y0 + m0 ( x − x0 ) x ∈[ x0 , x1 ]
 y + m (x − x ) x ∈[ x1 , x2 ]
 1 1 1
S ( x) = 
 ... ...

 yn −1 + mn −1 ( x − xn −1 ) x ∈[ xn −1 , xn ]

En la figura 2.24, se presenta el programa InteSpl1.m, para la

interpolación segmentaria; en la figura 2.25, una ejecución para los datos del
problema anterior.

f(x)
(x3,y3)

(x1,y1)

(x0,y0) (x2,y2)

x0 x1 x2 x3 x
Figura 2.23

function p=InteSpl1(x,y,x0)
% Interpolacion segmentaria spline lineal
% E. Raffo Lecca
%
% Datos
% n = numero de puntos
% x = vector x
% y = vector y
% x0 = valor a interpolar
% Resultados
 95

% p = valor interpolado
n = length(x);
if x0<x(1)
p=y(1);
else
if x0<=x(n)
i=1;
while x0>x(i)
i=i+1;
end
if x0==x(i)
p=y(i);
else
m=(y(i)-y(i-1))/(x(i)-x(i-1));
p=y(i-1)+m*(x0-x(i-1));
end
else
p=y(n);
end
end

Figura 2.24: InteSpl1.m

2.4.2 Spline cúbico

Sea el intervalo [a,b] con (n+1) puntos con:
a = x0 < x1 < x2 <...< xn = b
96

Figura 2.25: Ejecución de InteSpl1.m

La función S es llamada un spline cúbico, si cumple:

1. Es doblemente diferenciable y continua en [a,b].

2. S coincide en cada intervalo [x i, x i+1], i = 0,1,...,(n-1), con un
polinomio de grado3.

La función S spline cúbico posee n polinomios cúbicos S(x), con las

siguientes propiedades:

S i ( x ) = αi + βi ( x − xi ) + γ i ( x − xi ) 2 + δi ( x − xi ) 3
1.
para : x ∈[ xi , xi +1 ] , i = 0, 1, 2,..., n

2. S(x i) = y i el spline pasa por cada punto, i = 0,1,..., n.

3. S i(x i+1) = S i+1 (x i+1) el spline forma una función continua, i =
0, 1,..., n-2
4. S´i (x i+1) = S´i+1(x i+1) el spline forma una función suavisada, i
= 0, 1,..., n-2
 97

5. S´´ i (x i+1) = S´´ i+1(x i+1) la segunda derivada es continua, i = 0,

1,..., n-2

Para cada partición en [a,b], se tiene :

hi+1 = xi+1 – xi i = 0,1,..., n-1

Mi = S´´(xi ), siendo Mi el momento de S.

De la definición de la interpolación de Lagrange y sabiendo que la derivada

segunda en una función cúbica, coincide con una función lineal en cada
intervalo; se tiene:

x − xi +1 x − xi
S i´´ ( x ) = S ´´ ( xi ) + S ´´ ( xi +1 )
xi − xi +1 xi +1 − xi
x − xi +1 x − xi
S i´´ ( x ) = M i + M i +1 para x ∈[ xi , xi +1 ]
xi − xi +1 xi +1 − xi

al integrar S´´i (x), se tiene:

( xi +1 − x ) 2 ( x − xi ) 2
S i´ ( x ) = −M i + M i +1 + Ai
2 hi +1 2 hi +1
( x − xi +1 ) 3 ( x − xi ) 3
S i ( x) = M i + M i +1 + Ai ( x − xi ) + Bi
6 hi +1 6 hi +1
con i = 0, 1, 2, ..., n − 1

Desde Si (xi ) = yi y Si (xi+1) = yi+1, se tiene

2
h
M i i +1 + Bi = yi
6
2
h
M i +1 i +1 + Ai hi +1 + Bi = yi +1
6
98

Las constantes de integración Ai, Bi son:

2
h y − y i hi +1
Bi = y i − M i i +1 , Ai = i +1 − ( M i +1 − M i )
6 hi +1 6

Para determinar los coeficientes α i , βi , γ i y δ , se plantea:

S i ( xi ) = y i = α i , α i = yi
Mi
S i´´ ( xi ) = 2γ i = M i , γi =
2
M i hi +1 y − y i 2 M i + M i +1
S i´ ( xi ) = − + Ai = β i , β i = i +1 − hi +1
2 hi +1 6
M i +1 − M i M i +1 − M i
S i´´ ( xi ) = 6 δ i = , δi =
hi +1 6hi +1

Por cada [xi, xi+1] existen cuatro incógnitas, luego el spline cúbico
tiene aparentemente 4n incógnitas. Se plantean las siguientes ecuaciones:

De la propiedad 4, S’i (xi+1 ) = S’i+1(xi+1), con :

 99

− M i ( xi +1 − x) 2 ( x − xi ) 2 y i +1 − y i
S i´ ( x) = + M i +1 + −
2 hi +1 2 hi +1 hi +1
hi +1
( M i +1 − M i )
6
M h y − y i M i hi +1
S i´ ( xi +1 ) = i +1 i +1 + i +1 +
3 hi +1 6
− M i +1 ( x − xi +2 ) 2 M i +2 ( x − xi +1 ) 2 y i +2 − y i +1
S i´+1 ( x) = + + −
2 hi +2 2 hi +2 hi +2
hi +2
( M i +2 − M i +1 )
6
M h y − y i +1 M i +2 hi +2
S i´+1 ( xi +1 ) = − i +1 i +2 + i +2 −
3 hi +2 6

Se tiene un conjunto de ecuaciones:

hi M i −1 hi + hi +1 h y − y y − yi −1
+ M i + i +1 M i +1 = i +1 i − i
6 3 6 hi +1 hi
i = 1,2,..., n − 1

Como α, β, γ, δ , son función de M0, M1,..., Mn, entonces sólo se requieren

(n+1) incógnitas. Se tienen (n-1) ecuaciones y las otras dos ecuaciones se
encuentran cuando se cumple:

S0´´(a) = M0 = 0 = Mn = Sn-1´´(b)

en este caso se dice que el spline es natural, cuando:

S0´(a) = y0´ y Sn-1´(b) = yn´ ,

y se tiene las ecuaciones:

100

h1 h y −y
M 0 + 1 M 1 = 1 0 − y0´
3 6 h1
hn h y − yn −1
M n −1 + n M n = y´n − n
6 3 hn

en el caso del spline tipo grampa ó empalme (clamp).

Introduciendo las siguientes abreviaturas:

hi +1
λi = , μ i = 1 − λi
hi + hi +1

6  yi +1 − yi yi − yi −1 
bi =  −  i = 1, 2,..., n − 1
hi + hi +1  hi +1 hi 

con:
6  y1 − y 0 
λ 0 = 1 , b0 =  − y 0´ 
h1  h1 
6 y − y n −1 
μ n = 1 , bn =  y ´n − n 
hn  hn 

para el spline empalme. λ 0 = 0 , b0 = 0 , µ n = 0 , bn = 0 es para el spline

natural. Se tiene el sistema de ecuaciones lineales:

2M 0 + λ 0 M 1 = b0
μ1M 0 + 2 M 1 + λ1M 2 = b1
...
μ n −1M n −2 + 2 M n −1 + λ n −1M n = bn −1
μ n M n −1 + 2M n = bn

En notación matricial:
 101

2 λ0 0 ... 0 0 0   M 0   b0 
    
 μ1 2 λ1 ... 0 0 0   M 1   b1 
 ... ... ... ... ... ... ...   ...  =  ... 
    
0 0 0 ... μ n −1 2 λ n −1   M n −1   bn −1 
    
0 0 0 ... 0 μn 2   M n   bn 

Que es un sistema tridiagonal (ver capítulo 4) con:

ai = μ i
di = 2
ci = λ i , i = 0,1,2,..., n

Con a0 = 0 y cn = 0.

Para n =3 y los puntos (x,y) : (0,0), (1, 0.5), (2,3), (3,2) ; con la
modalidad natural se tienen los resultados siguientes:

x α β Γ Δ
(x0, x1) 0 -0.267 0 0.767
(x1, x2) 0.5 2.033 2.300 -1.833
(x2, x3) 3.0 1.133 -3.200 1.067

function p=InteSpl3(X,Y)
% Interpolacion segmentaria (spline) lineal
% E. Raffo Lecca
%
% Datos
% n =numero de puntos
% X =vector x
% Y =vector y
% x0=valor a interpolar
% Resultados
% p =valor interpolado
n=length(X);
while 1
option=input('Ingres opci¢n(1=Nor,2=Empal) -> ');
if (option>=1 & option<=2)
102

break;
end
end
if option==2
Y0Deriv=input('ingrese YDeriv[0] -> ');
YNDeriv=input('ingrese YDeriv[N] -> ');
end

%SplineCubic
for i=1:n-1
H(i+1)=X(i+1)-X(i);
end
for i=2:n-1
C(i)=H(i+1)/(H(i)+H(i+1));
D(i)=2;
A(i)=1-C(i);
b(i)=6*((Y(i+1)-Y(i))/H(i+1)-(Y(i)-Y(i-1))/H(i))/
(H(i)+H(i+1));
end
A(1)=0;
D(1)=2;
C(n)=0;
D(n)=2;
if option==1
b(1)=0;
C(1)=0;
b(n)=0;
A(n)=0;
else
C(1)=1;
b(1)=6*((Y(2)-Y(1))/H(2)-Y0Deriv)/H(2);
A(n)=1;
b(n)=6*(YNDeriv-(Y(n)-Y(n-1))/H(n))/H(n);
end

%TriDiagonal
M=Tridiagonal(A,D,C,b);
for i=1:n-1
Alpha(i)=Y(i);
Betha(i)=(Y(i+1)-Y(i))/H(i+1)-(2*M(i)+M(i+1))
*H(i+1)/6;
Gamma(i)=M(i)/2;
Delta(i)=(M(i+1)-M(i))/(6*H(i+1));
end
 103

% EscribirVector
fprintf('\n');
fprintf(' Alpha Betha Gamma
Delta\n');
for i=1:n-1
fprintf('%10.3f%10.3f%10.3f%10.3f\n',Alpha(i),Betha(i),...
Gamma(i),Delta(i));
end
end

Figura 2.26: InteSpl3.m

Figura 2.27: Ejecución de InteSpl3.m

En la figura 2.26, se presenta el programa InteSpl3.m, para la

interpolación segmentaria cúbica; en la figura 2.27, la ejecución para los datos
del problema anterior.
104

PROBLEMAS
1. Asumir f continua en [a,b] y x0 es un valor en Î [a,b]

f (x) = Pn(x) + Rn(x)

Donde:
∞
f ( k ) ( x0 )
Pn ( x) = ∑ ( x − x0 ) k
k =0 k!
f ( n +1)
Rn ( x) = ( x − x0 ) n +1
(n + 1) !
x3 x5 x7 x9
a) Si |x| < 1 la aproximación: sen ( x) = x − + − +
3! 5! 7! 9!

Encontrar P5(0.8) , P7(0.8) , R5(0.8) , R7(0.8)

b) Plantear un programa que evalúe sen(x).

2. Definiendo
n
P ( x) = ∏ ( x − xi )
i =0
Probar que:

Lk ( x) = ∏ ( x)
( x − x )∏ ´(x )
k k

3. Usar los polinomios de Lagrange para calcular

y(1.18), desde la tabla siguiente:
x y
1.05 1.025
1.10 1.049
1.15 1.073
1.20 1.095
1.25 1.118
1.30 1.140
 105

4. Usar la interpolación de Newton para calcular y(15),

desde la tabla :

i xi yi
0 6 104
1 10 160
2 20 370

5. Una serie de la forma:

Tm ( x ) =
a0 m
[
+ ∑ a j cos( jx ) + b j sen ( jx ) ,
2 j =1
]

es llamada una aproximación polinomial trigonométrica de orden n ;

siendo el grado m ≤ (n-1)/2.

Los coeficientes aj y bj son calculados por las fórmulas:

2 n
aj = ∑ f ( xi ) cos( jxi )
n i =1
, j = 0,1,..., m

2 n
bj = ∑ f ( xi ) sen( jxi )
n i =1
, j = 1,2,..., m

Asumir que los pares (xi , yi) están definidos como :

yi = f ( xi )
xi = − π + 2 iπ / n i = 1,2,..., n

6. Utilizar la función griddata de MATLAB, para generar variables aleatorias

(matrices x, y)
x ~ U(-2,2)
y ~ U(-2,2)
2
− y2
z = xe − x
106

7.
Utilizar la función Lagrange.m para ilustrar el uso de este método.

8. Demostrar la identidad Tm +n ( x ) +Tm −n ( x) = 2Tm ( x )Tn ( x) .

9. Presentar un M-file, para determinar las raíces de los polinomios de

Chebyshev; usando únicamente operadores matriciales de MATLAB.

10. Presentar un M-file, para determinar las raíces de los polinomios de

Chebyshev; usando la relación de recurrencia Tn(x).

11. Otra forma de presentar la interpolación, es conocer los coeficientes del

polinomio.

Sea el polinomio de interpolación escrito como:

y = a0 + a1 x + a2 x 2 + ... + an x n

Donde existen (n+1) pares de puntos (x,y); y se desea conocer los (n+1)
coeficientes de a. La solución al sistema de ecuaciones siguientes:

1 x0 x0 ... x0  a0   y0 
2 n
    
1 x1 x1 ... x1  a1   y1 
2 n

  =
 ... ... ... ...  ...   ... 
1 xn  an   yn 
2 n
 xn xn

Es conocida como matriz Vandermonde. Presentar el código MATLAB

para este problema.

Sugerencia:

En Numerical Recipies , se presenta el algoritmo de G. Rybicki; que utiliza

la derivada de un polinomio P(x), para encontrar los coeficientes.

12. En una interpolación con los pares (x, y), se quiere encontrar y(x) = c. El
caso de calcular x(y), y = c, éste es conocido como interpolación inversa.
 107

Los roles de x e y, son intercambiados. Un caso típico, es utilizarlo, para

computar la raíz de una función (en este caso y=0). Para la siguiente tabla,
determinar la raíz.

xk yk
1 -4
2 -2
4 2
6 1

13. Preparar el M-file, para el gráfico de la figura 2.11.

RESPUESTAS
1. P5(x) = x – x3 / 3! + x5 / 5! = 0.71739733

P5(x) = – x7 / 7! = - 0.00004161
P7(x) = x – x3 / 3! + x5 / 5! – x7 / 7! = 0.71735572
P7(x) = x9 / 9! = 0.0000003

n
2. De P ( x) = ∏ ( x − xi )
i =0

n
Luego P ´(x) = ∏ ( x − xi ) , y
i≠k

n
P ´(xk ) = ∏ ( xk − xi )
i≠k
De la definición del polinomio:

∏ (x − x ) i
L ( x) = i≠k

∏ (x − x )
k
k i
i≠k
Reemplazando:
108

∏ ( x − x )( x − x )
i k
∏ ( x)
L ( x) = i≠k
=
∏ ( x − x )( x − x ) ( x − x )∏´(x )
k
k i k k k
i≠k

3. Aplicando Lagrange, se tiene P(x) =1.08623.

i f (xi) Li(x)
0 1.024 0.01075
1 1.049 -0.08736
2 1.073 0.46592
3 1.095 0.69888
4 1.118 -0.09984
5 1.140 0.01165

4. Para las diferencias divididas finitas :

i xi F(xi ) 1 2
0 6 104 14 0.50
1 10 160 21
2 20 370

P(x) = 104 + 14(x-x0) + 0.50 (x-x0)(x-x1) = 252.50

5. Sea f(x) = x/2 sobre [-pi, pi], con n = 12 y m = 15. Los vectores resultantes:

a = (0.26180, -0.26180, 0.26180, -0.26180, 0.26180, -0.26180) , y

b = (0, 0.97705, -0.45345, 0.26180, -0.15115, 0.07015)

6.

function griddataDemo
% Demo para gridData
% E. Raffo Lecca
x =rand(100,1)*4-2;
y =rand(100,1)*4-2;
z =x.*exp(-x.^2-y.^2);
t =-2:0.20:2;
[XI,YI] =meshgrid(t,t);
ZI =griddata(x,y,z,XI,YI);
 109

mesh(XI,YI,ZI), hold
plot3(x,y,z,'o'), hold off

Figura 2.28: griddataDemo.m

Figura 2.29: Ejecución de griddataDemo.m

7.

function LagranDemo()
% Demo Interpolacion de Lagrange
% E. Raffo Lecca
x=[0 1 4 6]; y=[1 -1 1 -1];
x0=0:0.10:6;
for i=1:length(x0)
y0(i)=lagrange(x,y,x0(i));
end
plot(x,y,'o');
hold on;%para colocar varios plot;por default es off
110

plot(x0,y0);
legend('dato','Lagrange');
title('Interpolación de Lagrange');

Figura 2.30: lagranDemo.m

Figura 2.31: Ejecución de lagranDemo.h

8. Desde la identidad: 2 cos a cos b = cos( a + b) + cos( a − b) , se tiene:

2 cos mθ cos nθ = cos( m + n)θ + cos( m − n)θ , equivalente a

2Tm ( x )Tn ( x ) =Tm+n ( x ) +Tm−n ( x )

9.

function x=rootCheby1(n)
 111

% Raices de polinomio de Chebyshev

% E. Raffo Lecca
%
x=cos(0.5*(2*(0:n-1)+1)*pi/n);

Figura 2.32: M-file rootCheby1.m

Comentario:

2*(0:n-1)+1 , es equivalente a (2k-1), k = 1, ..., n.

10.

function x=rootCheby2(n)
% Raices de polinomio de Chebyshev
% E. Raffo Lecca
%
n=n+1;
if n==1
p=[1]; % p=1
else if n==2
p=[1,0];% p= 0+x
else
p2=[1];
p1=[1,0];
for k=3:n
p=2*[p1,0]-[0,0,p2];
p2=p1;
p1=p;
end
end
end
x=sort(roots(p));

Figura 2.33: M-file rootCheby2.m

12.
P (0) = 5.67

k xk bk ck
3 2 -0.167 -0.167
2 1 0.167 0.344
 112

1 -2 0.5 1.168
0 -4 1 5.67

Figura 2.34: Ejecución para interploación inversa

13.
function polyCheby()
% E. Raffo Lecca
% Compara función soporte
% entre puntos equidistantes y chebyshev
n=10;
a=-1;b=1;
x=linspace(a,b,n+1);
xCheby=sort(cos(0.5*(2*(1:n)-1)*pi/n));
x0=a:0.01:b;
for k=1:length(x0)
w(k)=prod(x0(k)-x);
wCheby(k)=prod(x0(k)-xCheby);
end

plot(x0,w,x0,wCheby)
legend('equidistante','Chebyshev')

REFERENCIAS
1. F.B. Hildebrandt, Introduction to Numerical Analysis, McGraw-Hill Book
Company, Inc, 1956.
 113

2. G. Dahlquist, A. Bjorck, Numerical Methods in Scientific Computing

3. Williams Press , Saul A. Teukolsky, William Vetterling, Brian P.
Flannery, Numerical Recipies in C++, The Art of Scientific
Computing, Cambridge University Press, 2002, Second Edition.
4. Eduardo Raffo Lecca, MATLAB: Mi Primer, Raffo-Lecca editors
2001.
5. J. Stoer, R. Bulirsch., Introduction to Numerical Analysis, Springer
-Verlag, 1993, Second Edition.
6. Wong, Samuel S. M., Computational Methods in Physics and
Engineering, World Scientific Publishing Co. Pte. Ltd, 1997, Second
Edition.