Interpolación y aproximación polinomial
Programación Numérica
Definición
Un polinomio de grado n es una expresión de la forma:
P(x) = anxn + an-1xn-1 + ... +a1x + a0
Donde an <> 0
Teorema (teorema de aproximación de Weierstrass)
Suponga que f está definida y es continua en [a, b]. Para > 0 existe un polinomio P definido en [a, b], con la propiedad de que
|f(x) – P(x)| < , para toda x en [a, b]
Desarrollo en series de TaylorSea f(x) = ex
Desarrollando en serie de Taylor alrededor de x = 0
P0(x) = 1 P1(x) = 1 + x P2(x) = 1 + x + x2/2
P3(x) = 1 + x + x2/2 + x3/6 P4(x) = 1 + x + x2/2 + x3/6 + x4/24
P5(x) = 1 + x + x2/2 + x3/6 + x4/24 + x5/120
Valores de ex
x p0(x) p1(x) p2(x) p3(x) p4(x) p5(x) exp(x)-2.0 1.00000 -1.00000 1.00000 -0.33333 0.33333 0.06667 0.13534-1.5 1.00000 -0.50000 0.62500 0.06250 0.27344 0.21016 0.22313-1.0 1.00000 0.00000 0.50000 0.33333 0.37500 0.36667 0.36788-0.5 1.00000 0.50000 0.62500 0.60417 0.60677 0.60651 0.606530.0 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.000000.5 1.00000 1.50000 1.62500 1.64583 1.64844 1.64870 1.648721.0 1.00000 2.00000 2.50000 2.66667 2.70833 2.71667 2.718281.5 1.00000 2.50000 3.62500 4.18750 4.39844 4.46172 4.481692.0 1.00000 3.00000 5.00000 6.33333 7.00000 7.26667 7.38906
Valores de las aproximaciones de ex con polinomios de Taylor
Expansión de Taylor para 1/x
n
k
kkn
k
kk
n xxk
fxP
00
111!1
n 0 1 2 3 4 5 6 7
Pn(3) 1 -1 3 -5 11 -21 43 -85
Interpolación polinomial de Newton
Revisaremos solo algunos casos: lineal, de segundo grado y de tercer grado.
Interpolación lineal
x0 x x1
f(x0)
f1(x)
f(x1)
f(x)
Utilizando triángulos semejantes
01
01
0
01
xxxfxf
xxxfxf
Reordenando
001
0101 xx
xxxfxf
xfxf
EjemploEstimar ln 2 mediante interpolación lineal si ln1 = 0 y ln 6 = 1.791759 y ln 4 = 1.386294
001
0101 xx
xxxfxf
xfxf 3583519.012
160791759.1
1ln21
f
4620981.01214
0386294.11ln21
f
Valor real ln 2 = 0.6931472
Error relativo porcentual = 33.3%
0 1 2 3 4 5 6 7 8-1
-0.5
0
0.5
1
1.5
2
2.5 f(x) = ln x
f1(x)
Estimaciones lineales
Valor verdadero
Interpolación cuadráticaPolinomio cuadrático
f2(x) = b0 + b1(x – x0) + b2(x – x0)(x – x1) (1)
simplificado
f2(x) = b0 + b1x – b1x0 + b2x2 + b2x0 x1 – b2xx0 – b2xx1
Podemos escribirlo como
f2(x) = a0 + a1x + a2x2
Donde
a0 = b0 – b1x0 + b2x0 x1, a1 = b1 – b2x0 – b2x1, a2=b2
Podemos evaluar b0, b1 y b2 sustituyendo x0, x1 y x2 en la ecuación (1), se obtiene
b0 = f(x0)
01
011 xx
xfxfb
02
01
01
12
12
2 xxxx
xfxfxx
xfxf
b
ejemplo 2Calculemos ln 2 con ln 4 y ln 6, los punto que se conocen son:
x0 = 1 f(x0) = 0
x1 = 4 f(x0) = 1.386294
x0 = 6 f(x0) = 1.791759
Aplicando las ecs. anteriores
b0 = 0
b1 = (1.386294 – 0)/(4 – 1) = 0.4620981
b2 = ((1.791759 – 1.386294)
/(6 – 4) – 0.4620981)/(6 – 1)
= – 0.0518731
El polinomio es
f2(x) = 0.4620981(x – 1) – 0.0518731(x – 1)(x – 4)
f2(2) = 0.5658444
0 1 2 3 4 5 6 7 8-1
-0.5
0
0.5
1
1.5
2
2.5
f(x) = ln x
Estimación cuadrática
Valor verdadero
Estimación lineal
Valor real ln 2 = 0.6931472
Error relativo porcentual = 18.4%
Forma generalPolinomio general
fn(x) = b0 + b1(x – x0) +...+ bn(x – x0)(x – x1)... (x – xn–1)
Los coeficientes se calculan con
b0 = f(x0)
b1 = f [x1, x0]
b2 = f [x2, x1, x0]
bn = f [,xn, xn–1, ..., x1, x0]
Donde los paréntesis cuadrados se denominan diferencias divididas finitas.
La n-ésima diferencia dividida finita es:
Se conoce como polinomio de interpolación de Newton en diferencias divididas.
0
02111011
,...,,,...,,,,...,,
xxxxxfxxxf
xxxxfn
nnnnnn
ejemplo 3Calculemos ln 2 con ln 0, ln 4, ln 5 y ln 6, los punto que se conocen son:
x0 = 1 f(x0) = 0
x1 = 4 f(x1) = 1.386294
x2 = 6 f(x3) = 1.791759
x3 = 5 f(x2) = 1.609438
primeras diferencias
f [x1, x0] = (1.386294 – 0)/(4 – 1) = 0.4602981
f [x2, x1] = (1.791759 – 1.386294)/(6 – 4) = 0.2027326
f [x3, x2] = (1.609438 – 1.791759)/(5 – 6) = 0.1823216
Segundas diferencias
f [x2, x1, x0] = (0.2027326 – 0.4602981)/(6 – 1) = –0.05187311
f [x3, x2, x1] = (0.1823216 – 0.2027326)/(5 – 4) = –0.02041100
tercera diferencia
f [x3, x2, x1 , x0] = (–0.02041100–(–0.05187311))/(5 – 1) = 0.007865529
Polinomio
f3(x) = 0 + 0.4602981(x – 1) –0.05187311(x – 1) (x – 4) + 0.007865529(x – 1) (x – 4) (x – 6)
Valor calculado con el polinomio
f3(2) = 0.6287686
Ejemplo 3 (cont.)
0 1 2 3 4 5 6 7 8-1
-0.5
0
0.5
1
1.5
2
2.5
f(x) = ln x
Valor verdadero
Estimación cúbica
f3(x)
Estimación del error
Rn = f [,xn+1, xn, ..., x1, x0](x – x0) (x – x1)... (x – xn)
Para estimar el error requerimos de un datos más (xn+1). La siguiente fórmula puede utilizarse para estimar el error.
Interpolación y polinomio de Lagrange
Se trata de encontrar un polinomio de grado n que pase por los puntos (x0, f(x0)), (x1, f(x1)), ... (xn, f(xn)), se construye un cociente Ln,k(xk) con la propiedad de que
Ln,k(xi) = 0 cuando i k y Ln,k(xk) = 1
Se requiere entonces que el numerador contenga
(x – x0) (x – x1)... (x – xk–1)(x – xk+1)... (x – xn)
El denominador debe coincidir con el numerador cuando x = xk.
n
kii ik
i
nkkkkkkk
nkkkn xx
xxxxxxxxxxxx
xxxxxxxxxxxL
01110
1110,
N-ésimo polinomio interpolante de Lagrange
Teorema
Si x0, x1, x2, ... xn, son n+1 números distintos y si f es una función cuyos valores están dados en esos números, entonces existe un polinomio de grado a lo más n, con la propiedad de que
f(xk) = P(xk) para cada k = 0, 1, 2, ...n
Este polinomio está dado por
n
kknknnnn xLxfxLxfxLxfxP
0,,0,0
n
kii ik
i
nkkkkkkk
nkkkn xx
xxxxxxxxxxxx
xxxxxxxxxxxL
01110
1110,
donde
Aproximación a 1/x con interpolantes de Lagrange
P(x) = 0.5*((x–6.5)x+10)+0.4*((–4x+24)x–32)/3+ 0.25*((x + 4.5)x+5)/3
P(x) = (0.05x – 0.425)x + 1.15 = 0.05x2 – 0.425x + 1.15
f(3) = P(3) = 0.325
Usaremos x0 = 2, x1 = 2.5 y x2 = 4, para obtener un polinomio de grado 2 para 1/x. f(x0) = 0.5, f(x1)= 0.4 y f(x2) = 0.25.
Los polinomios de Lagrange son:
3
55.45.24245.22
2,
332244
45.225.242
1,
105.6425.0245.2
0,
xxxxx
nL
xxxxx
nL
xxxx
xn
L
Aproximación a 1/x con interpolantes de Lagrange
P(x) = (0.05x – 0.425)x + 1.15
f(3) = P(3) = 0.325
El error en la interpolación de Lagrange
El error en la interpolación de Lagrange puede calcularse con
n
n
xxxxxxn
xfxPxf
...!1 10
1
0
Algoritmo en Matlabfunction fi = Lagran_(x, f, xi)
fi=zeros(size(xi));
np1=length(f);
for i=1:np1
z=ones(size(xi));
for j=1:np1
if i~=j, z = z.*(xi - x(j))/(x(i)-x(j));end
end
fi=fi+z*f(i);
end
return
Calcula coeficientes de P2(x)
%Calcula el polinomio interpolante de Lagrange de grado 2
function [a,b,c] = Lagrange(x0,x1,x2,fx0,fx1,fx2)
t0 = (x0 - x1)*(x0 - x2);
t1 = (x1 - x0)*(x1 - x2);
t2 = (x2 - x0)*(x2 - x1);
a = fx0/t0 +fx1/t1 +fx2/t2;
b = -fx0*(x1 + x2)/t0 - fx1*(x0 + x2)/t1 - fx2*(x0 + x1)/t2;
c = fx0*x1*x2/t0 + fx1*x0*x2/t1 + fx2*x0*x1/t2;
Interpolación InversaTabla de valores de f (x) = 1/x.
x 1 2 3 4 5 6 7
f (x) 1 0.5 0.3333 0.25 0.2 0.1667 0.1429
Se desea conocer el valor de x tal que f (x) = 0.3.
El problema se resuelve definiendo un polinomio de interpolación de grado 2 con los puntos (2, 0.5), (3, 0.3333) y (4, 0.25) y resolviendo la ecuación:
f (x) = 0.3 = 1.08333 – 0.375x + 0.041667x2
Lo que da x = 5.704158 y x = 3.295842, el valor real es 3.333.
Trazadores (Splines)Dados n +1 puntos podemos construir un polinomio de grado n para interpolar valores dentro del intervalo.
También se pueden usar líneas rectas entre cada par de puntos para hacer interpolación lineal entre ellos o polinomios cuadráticos o cúbicos.
Tales interpoladores se llaman trazadores lineales, cuadráticos y cúbicos, respectivamente.
La ventaja de los trazadores es que no presentan el efecto de oscilación de los polinomios de alto grado.
f (x)
x
f (x)
f (x)
f (x)
Trazadores linealesPara los trazadores lineales se definen rectas entre cada intervalo para calcular los valores intermedios.
f (x) = f (x0) + m0(x – x0) x0 <= x <= x1
f (x) = f (x1) + m1(x – x0) x1 <= x <= x2
f (x) = f (x0) + mn–1 (x – x0) xn–1 <= x <= xn
Los valores de mi se calculan con:
ii
iii xx
xfxfm
1
1
ejemplo
x f (x)
3.0 2.5
4.5 1.0
7.0 2.5
9.0 0.52 4 6 8 10
2
0
Trazadores cuadráticosEl polinomio en cada intervalo es de la forma:
fi(x) = ai x2 + bi x + ci
Para encontrar los ai , bi , ci se deben cumplir las siguientes condiciones:
1. Los valores de la función deben ser iguales en los nodos interiores, 2n – 2 ecuaciones.
2. La primera y última función debe pasar por los extremos, 2 ecuaciones.
3. Las primeras derivadas en los nodos interiores deben ser iguales, n – 1 ecuaciones. O sea: 2ai –1 xi–1 + bi –1 = 2ai xi–1 + bi
4. Suponer derivada 0 en el primer punto. a1 = 0
)(
)(
112
1
11112
11
iiiiii
iiiiii
xfcxbxa
xfcxbxa
)(
)(2
0101201
nnnnnn xfcxbxa
xfcxbxa
ejemplo
x f (x)
3.0 2.5
4.5 1.0
7.0 2.5
9.0 0.5
Encontrar f (5)
La condición 1 genera las siguientes ecuaciones:
20.25a1 + 4.5b1 + c1 = 1.0
20.25a2 + 4.5b2 + c2 = 1.0
49a2 + 7b2 + c2 = 2.5
49a3 + 7b3 + c3 = 2.5
La condición 2 da las siguientes ecuaciones
9a1 + 3b1 + c1 = 2.5
81a3 + 9b3 + c3 = 0.5
La condición 3 genera:
9a1 + b1 = 9a2 + b2
14a2 + b2 = 14a3 + b3
0
0
5.0
5.2
5.2
5.2
1
1
0114011400
00001901
198100000
00000013
174900000
000174900
00015.425.2000
00000015.4
3
3
3
2
2
2
1
1
c
b
a
c
b
a
c
b
El sistema resultante es:
La solución es:
a1 = 0 b1 = –1 c1 = 5.5
a2 = 0.64 b2 = –6.67 c2 = 18.46
a3 = –1.6 b3 = –24.6 c3 = 91.3
f(5) = 0.64(5)2 – 6.67(5) +18.46 = 1.11
yi = ppval (pp, xi) - Evalúa polinomio a trozos pp en los puntos xi. Si pp.d es un escalar mayor que 1, o un arreglo, entonces el valor regresado yi será un arreglo que es d1, d1, ..., dk, length (xi).
pp = spline (x, y) yi = spline (x, y, xi)Regresa los interpolantes cúbicos de y en los puntos x. Si se llama con dos argumentos, regresa los trozos polinomicos pp que ueden ser evaluados con ppval.Si se llama con tres parámetros, evalúa el los puntos xi.
Splines cúbicos
xxhz
hy
xxhz
hy
xxh
zxx
hz
xS
iii
i
ii
ii
i
i
ii
ii
i
ii
111
3131
66
66
Aplicando las condiciones de continuidad del spline S y de las derivadas primera S' y segunda S'', es posible encontrar la expresión analítica del spline.
Donde hi = xi+1 – xi y z0, z1, … ,zn son incognitas.
Aplicando las condiciones de continuidad se llega a
11
11
1111
662
ii
iii
iiiiiiii yy
hyy
hzhzhhzh
La ecuación anterior, genera un sistema de n–1 ecuaciones lineales con n+1 incógnitas.
2
3
2
1
2
3
2
1
12
22
32
221
11
0000
000
000
000
000
0000
n
n
n
n
n
n
v
v
v
v
z
z
z
z
uh
hu
uh
huh
hu
1
111
1
1
21
1
6
2
i
iiiii
iii
i
i
iiii
u
vhbbv
yyh
b
u
hhhu
donde
Los valores del spline S se calculan eficientemente con
iiiiiiii AxxBxxCxxyxS
Donde
iii
ii
ii
i
ii
iii
i
yyh
zh
zh
C
zB
zzh
A
11
1
136
2
61
Los coeficientes de los polinomios se pueden calcular con:c1 = yi – xi Dc2 = D – xi Ec3 = E – xi Ac4 = A
Para obtener: fi (x) = c1 + c2 x + c3 x2 + c4 x3
Donde
AxBE
AxBxCD
yyh
zh
zh
C
zB
hzz
A
i
ii
iii
ii
ii
i
i
ii
2
136
2
6
2
11
1
Guión en MatLab%encuentra los trazadores cúbicos para un conjunto de puntos x,y% x - vector con los n valores de x% y - vector con los n valores de y% w - matriz de n-1 por 4 con los coeficientes de los polinomios cúbicosfunction w = spline3(x,y) [dummy n] = size(x); for i = 1:n-1 h(i) = x(i+1)-x(i); b(i) = 6*(y(i+1)-y(i))/h(i); end u(2) = 2*(h(1)+h(2)); v(2) = b(2)-b(1); for i = 3:n-1 u(i) = 2*(h(i)+h(i-1))-h(i-1)^2/u(i-1); v(i) = b(i)-b(i-1)-h(i-1)*v(i-1)/u(i-1); end
z(n) = 0; for i = n-1:-1:2; z(i) = (v(i)-h(i)*z(i+1))/u(i); end z(1) = 0; for i = 1:n-1 A = (z(i+1)-z(i))/6/h(i); B = z(i)/2; C = -h(i)*z(i+1)/6-h(i)*z(i)/3+(y(i+1)-y(i))/h(i); D = C-x(i)*B+A*x(i)^2; E = B-2*x(i)*A; w(i,4) = y(i)-x(i)*D; w(i,3) = D-x(i)*E; w(i,2) = E-x(i)*A; w(i,1) = A; endend
Ejemplo