Diferencias divididas e interpolación osculatoria

Traballo Fin de Grao Diferencias divididas e interpolación osculatoria Elba García Hermida 2018/2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Diferencias divididas e interpolación osculatoria Elba García Hermida Julio de 2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Matemática Aplicada Título: Diferencias divididas e interpolación osculatoria Breve descrición do contido El trabajo contará con una parte teórica y una práctica. En la parte teórica se abordarán los contenidos siguientes: (i) Diferencias divididas: definición, fórmula recursiva en el caso de que todos los argumentos sean distintos y expresión del polinomio de interpolación de Lagrange mediante diferencias divididas (fórmula de interpolación de Newton). (ii) Interpolación osculatoria (interpolación de Hermite generalizada): concepto, existencia y unicidad del polinomio de interpolación osculador, expresión del error. (iii) Diferencias divididas con argumentos repetidos. Expresión del polinomio de interpolación osculador mediante diferencias divididas. La parte práctica consistirá en programar en matlab algunos métodos de cálculo del polinomio de interpolación osculador descritos en la parte téorica y aplicarlos a ejemplos concretos. iii Índice general Resumen vii Introducción ix 1. Diferencias divididas: Fórmula de interpolación de Newton 1 1.1. El polinomio de interpolación de Lagrange . . . . . . . . . . . . . . . . . . . 1 1.2. El polinomio de interpolación de Newton . . . . . . . . . . . . . . . . . . . . 5 1.2.1. Diferencias finitas: fórmulas de Newton progresiva y regresiva . . . . 10 2. Interpolación de Hermite 19 2.1. El polinomio de interpolación de Hermite . . . . . . . . . . . . . . . . . . . 19 2.2. Acotación del error para la interpolación de Hermite . . . . . . . . . . . . . 23 3. Diferencias divididas con argumentos repetidos 29 3.1. Diferencias divididas con argumentos repetidos . . . . . . . . . . . . . . . . 29 ANEXO I. Codigo Matlab para el polinomio de interpolación de Lagrange 39 ANEXO II. Codigo Matlab para el polinomio de interpolación de Hermite 43 Bibliografía 47 v Resumen En el primer capítulo demostraremos la existencia y unicidad del polinomio de interpolación de Lagrange y obtendremos la fórmula de Newton. Además, justificaremos el uso de las tablas de diferencias divididas para obtener la fórmula de Newton para el polinomio de interpolación. En segundo y tercer capítulo trataremos la existencia y unicidad del polinomio de interpolación de Hermite. Demostraremos también algunos resultados relativos al error para este tipo de interpolación. Extenderemos la definición de diferencia dividida al caso de argumentos repetidos y veremos que entonces el polinomio de Hermite se puede escribir en la forma de Newton. Abstract In the first chapter we will show the existence and uniqueness of the Lagrange’s interpolation polynomial. We will obtain the Newton’s formula. Furthermore, we will justify the use of the divided differences tables to find the Newton’s formula for the interpolation polynomial. In the second and third chapter we will deal with the existence and uniqueness of the Hermite’s interpolation polynomial, as well as some error related results. We will extend the divided difference definition to the case of repeated arguments. We will prove that we can write the Hermite’s polynomial in the Newton’s form. vii Definición 1.4. Sea f(x)una función dada y Pn(x)una aproximación polinómica de ésta; se define el error puntual entre f(x)yPn(x)por: Rn(x)≡f(x)−Pn(x)(1.7) Es de gran utilidad tener una expresión explícita para este error. El siguiente teorema nos proporciona una representación de Rn(x): Teorema 1.5. Supongamos que f(x)tiene derivada (n+1)-ésima, f(n+1)(x), en un intervalo [a, b]2. Sea Pn(x)el polinomio de interpolación de f(x)con respecto a n+ 1 puntos distintos xi,i= 0,1, ..., n, pertenecientes al intervalo [a, b]. Entonces, para cada x∈[a, b], x6=xi, i = 0,1, . . . , n, existe un punto ξ=ξ(x)en el intervalo abierto m´ın{x0, x1, ..., xn, x}< ξ < m´ax{x0, x1, ..., xn, x}(1.8) tal que f(x)−Pn(x)≡Rn(x) = (x−x0)(x−x1)...(x−xn) (n+ 1)! f(n+1)(ξ) ≡wn(x) (n+ 1)!f(n+1)(ξ) (1.9) Demostración: Definimos Sn(x), para cualquier x6=xi,i= 0, . . . , n, por: Rn(x)≡(x−x0)(x−x1)...(x−xn)Sn(x) = wn(x)Sn(x)(1.10) Consideramos x∈[a, b],x6=xi,i= 0, . . . , n, fijo y definimos la función F(z)por: F(z)≡f(z)−Pn(z)−wn(z)Sn(x).(1.11) Dado que ftiene derivada (n+ 1)-ésima en [a, b],F(n+1)(z)está definida en [a, b]. F(z)se anula en n+ 2 puntos distintos de [a, b]: F(x0) = F(x1) = ... =F(xn) = F(x) = 0. Por lo tanto, existen n+ 1 intervalos consecutivos contenidos en [a, b]en cuyos extremos F(z)es cero. En virtud del teorema de Rolle, en el interior de cada uno de estos intervalos existe al menos un punto en el cual F0(z)se anula. Tenemos, por tanto, al menos n+ 1 puntos distintos en el intervalo (1.8) en los cuales F0(z) = 0. A su vez, estos puntos definen al menos nintervalos en cuyo interior, aplicando de nuevo el teorema de Rolle, existe un punto en el que la derivada de F0(z)se anula. Con lo cual, F00(z)=0en al menos npuntos 2En realidad, la conclusión del teorema se obtiene del mismo modo para condiciones más débiles: basta que f(x)sea continua en el intervalo cerrado [a, b]y que tenga derivada (n+ 1)-ésima en el intervalo (a, b). 4 distintos en el intervalo (1.8). Repitiendo este proceso, llegamos a que existe un punto en el intervalo (1.8), que denotamos ξ, tal que F(n+1)(ξ) = 0. Puesto que Pn(z)es un polinomio de grado n, dn+1Pn dzn+1 (z)=0. Por otra parte, dn+1 dzn+1 [wn(z)Sn(x)] = (n+ 1)!Sn(x). Derivando en la expresión (1.11) y teniendo en cuenta las dos últimas ecuaciones se obtiene que: F(n+1)(z) = f(n+1)(z)−(n+ 1)!Sn(x) y puesto que F(n+1)(ξ)=0se deduce que: Sn(x) = 1 (n+ 1)!f(n+1)(ξ).(1.12) Sustituyendo (1.12) en (1.10) se concluye (1.9). Observación 1.6.Si x=xi,i= 0,1, . . . , n,Rn(x) = 0 y (1.9) también se verifica en este caso (cualquiera que sea ξ∈[a, b]). Si podemos acotar |f(n+1)(x)|en [a, b], mediante (1.9) obtendremos una acotación del error, como veremos en el Ejemplo 1.17 más adelante. 1.2. El polinomio de interpolación de Newton Esta sección recoge la obtención de la fórmula de Newton para el polinomio de interpolación de Lagrange, así como la introducción de las tablas de diferencias divididas. Examinaremos también el error cuando tomamos puntos uniformemente espaciados. Ya hemos visto que el polinomio de interpolación existe y es único; además, dado un conjunto de puntos de interpolación, lo podemos construir fácilmente usando los polinomios fundamentales de Lagrange. La principal desventaja de este método es que si añadimos un nuevo punto de interpolación, el nuevo polinomio de grado superior no puede ser obtenido simplemente modificando el anterior. Buscaremos una nueva representación del polinomio de interpolación que cuente con esta propiedad, es decir, que nos permita encontrar polinomios de interpolación de grado mayor de forma sencilla añadiendo nuevos términos. 5 Sea Pk(x)el polinomio de interpolación de f(x), de grado menor o igual que k, con respecto a los k+ 1 puntos distintos x0, x1, . . . , xk. Buscamos los polinomios de interpolación sucesivos {Pk(x)}, de grado menor o igual que ktales que P0(x)≡f(x0)y Pk(x) = Pk−1(x) + qk(x), para k = 1,2, . . . , n, (1.13) donde qk(x)tiene a lo sumo grado k. Puesto que necesitamos que Pk(xj) = f(xj) = Pk−1(xj), j = 0,1, . . . , k −1, deducimos que qk(xj)=0en estos kpuntos. Por ello, podemos escribir qk(x) = ak k−1 Y j=0 (x−xj),(1.14) expresión que representa el polinomio más general de grado menor o igual que kque se anula en estos kpuntos. Desconocemos el valor de las constantes akpero, de la ecuación Pk(xk) = f(xk), junto con (1.13) y (1.14), deducimos que ak=f(xk)−Pk−1(xk) Qk−1 j=0 (xk−xj), para k = 1,2, . . . , n. (1.15) Trivialmente, obtenemos el polinomio de interpolación de grado cero para el punto inicial x0:P0(x)≡f(x0). Definimos a0=f(x0). Usando repetidamente (1.13) y (1.14) para k= 1,2, . . . , n, obtenemos que Pn(x) = a0+a1(x−x0) + · · · +an(x−x0). . . (x−xn−1).(1.16) El coeficiente akrecibe el nombre de diferencia dividida de orden k; normalmente se emplea la notación: f[x0] := a0, f[x0, x1, ..., xk] := ak, k = 1,2, ... (1.17) Buscamos ahora una expresión alternativa para estos coeficientes, tratando que sea más explícita que la expresión (1.15). Dado que el polinomio Pn(x)obtenido en (1.16) es el único polinomio de interpolación de grado ≤n, tiene que ser el mismo que obtuvimos en forma de Lagrange en (1.5). Igualando el coeficiente principal en las expresiones (1.16) y (1.5), deducimos que: an=f[x0, x1, ..., xn] = n X j=0 f(xj) Qn k=0,k6=j(xj−xk).(1.18) 6 De la representación (1.18) se sigue que las diferencias divididas son funciones simétricas de sus argumentos. Esto es, si usamos la notación fi,j,k,... ≡f[xi, xj, xk, ...], entonces esta simetría se expresa por: f0,1,...,n =fj0,j1,...,jn,(1.19) siendo (j0, j1, ..., jn)una permutación arbitraria de los enteros (0,1, ..., n). Lema 1.7. f[x0, x1, . . . , xk] = f[x1, x2, . . . , xk]−f[x0, x1, . . . , xk−1] xk−x0 ,(1.20) con f[xi] = f(xi), i = 0,1,2, . . . , k. Demostración: Denotemos por P∗ k−1el polinomio de interpolación de fde grado ≤k−1 relativo a los puntos x1, x2, . . . , xk; el coeficiente de xk−1es f[x1, x2, . . . , xk]. Definimos el polinomio Qk(x) = ((x−x0)P∗ k−1(x)−(x−xk)Pk−1(x))/(xk−x0). Qktiene grado menor o igual que ky coincide con fen los puntos x0, x1, . . . , xk; por la unicidad del polinomio de interpolación, Qk(x)≡Pk(x). Igualando el coeficiente de xkde Qkcon el de Pkse obtiene (1.20)  Esta fórmula recursiva justifica el nombre de diferencias divididas y el uso de la tabla mostrada a continuación para obtener el polinomio de interpolación, pues cada diferencia dividida puede ser calculada explicitamente en función de f(x0), f(x1), . . . , f(xn) 7 Tabla 1. Diferencias Divididas. xf(x)f[x, x]f[x, x, x]f[x, x, x, x]. . . x0 x1 x2 x3 x4 . . . f0 f1 f2 f3 f4 . . . f1−f0 x1−x0≡f01 f2−f1 x2−x1≡f12 f3−f2 x3−x2≡f23 f4−f3 x4−x3≡f34 . . . f12−f01 x2−x0≡f012 f23−f12 x3−x1≡f123 f34−f23 x4−x2≡f234 . . . f123−f012 x3−x0≡f0123 f234−f123 x4−x1≡f1234 . . . . . . La primera columna está formada por los valores de la función f, la segunda por las diferencias divididas de primer orden, etc; pasamos de una columna a otra aplicando la fórmula (1.20). Así, el polinomio de interpolación se puede escribir como: Pn(x) = f[x0] + f[x0, x1](x−x0) + ... +f[x0, x1, ..., xn](x−x0)...(x−xn−1)(1.21) Esta expresión de Pn(x)se conoce como fórmula de interpolación de Newton; resulta muy sencilla de obtener a partir de la tabla. 8 Ejemplo 1.8. Calcularemos el polinomio de interpolación para la función ˆ K(ˆx)del Ejemplo 1.3 relativo a los puntos 2,3,5,6, en grados, empleando la Tabla 1. Obtenemos así la siguiente tabla de diferencias divididas: ˆxˆ K(ˆx)ˆ K[ˆx, ˆx]ˆ K[ˆx, ˆx, ˆx]ˆ K[ˆx, ˆx, ˆx, ˆx] 2 3 5 6 1,5713 1,5719 1,5738 1,5751 0,0006 0,00095 0,0013 1,1 Û 6·10−4 1,1 Û 6·10−4 0 La fórmula de Newton para el polinomio de interpolación de Lagrange de K(ˆx)relativo a los puntos 2,3,5,6 viene dada por: P3(ˆx)=1,5713 + 0,0006(ˆx−2) + 1,1 Û 6·10−4(ˆx−2)(ˆx−3) + 0(ˆx−2)(ˆx−3)(ˆx−5). P3(3,5) = 1,57228750 que es el valor que habíamos obtenido en el Ejemplo 1.3. Para realizar los cálculos hemos empleado el programa NEWTON.m, cuyo código se incluye en el anexo al final del documento. Nótese que, como ya se ha mencionado, para obtener un polinomio de interpolación de orden superior añadiendo un nuevo punto, tan solo necesitamos sumar otro término a (1.21) pero esta vez se ve involucrada una diferencia dividida de orden superior. De (1.16) y (1.17) se deduce que Pn+1(x) = Pn(x) + f[x0, . . . , xn, xn+1] n Y j=0 (x−xj) y como Pn+1(xn+1) = f(xn+1) f(xn+1)−Pn(xn+1) = f[x0, . . . , xn, xn+1] n Y j=0 (xn+1 −xj). 9 Sea x6=x0, . . . , xn, reemplazando xn+1 por xse obtiene f(x)−Pn(x) = f[x0, x1, ..., xn, x] n Y j=0 (x−xj),(1.22) Esta igualdad nos da una representación del error distinta de la vista en el Teorema 1.5. Proposición 1.9. Bajo las hipótesis del Teorema 1.5, para todo x∈[a, b], x 6=xi, i = 0,1, . . . , n, existe ξ=ξ(x)en el intervalo (1.8) tal que f[x0, x1, . . . , xn, x] = 1 (n+ 1)!f(n+1)(ξ).(1.23) Demostración: Basta comparar las representaciones del error (1.9) y (1.22).  1.2.1. Diferencias finitas: fórmulas de Newton progresiva y regresiva Cuando los puntos xide interpolación no son equidistantes, utilizamos la fórmula de Newton; sin embargo, cuando los xison equidistantes podemos emplear fórmulas más simples y menos costosas. Supongamos conocidos los valores de la función fen los puntos equidistantes xicon paso h > 0 xi=x0+ih, i = 0,1, . . . (1.24) Definición 1.10. Definimos el operador diferencia progresiva,4, por: 40fi=f(xi) = fi, 4fi=f(xi+1)−f(xi), y para todo entero k≥0, 4k+1fi=4(4kfi) = 4kfi+1 − 4kfi Lema 1.11. Se verifica que: f[xi, xi+1, . . . , xi+k] = 4kfi k!hk,∀k≥0. Demostración: Procederemos por inducción en k. Para k= 0 es trivial. Para k= 1, debido a (1.24), tenemos que f[xi, xi+1] = f(xi+1)−f(xi) xi+1 −xi =4fi h 10 Supongamos que el resultado se verifica para k=ny veamos que se cumple para k=n+ 1. Usando (1.20), la hipótesis de inducción y (1.24) tenemos que f[xi, xi+1, . . . , xi+n, xi+n+1] = f[xi+1, . . . , xi+n+1]−f[xi, . . . , xi+n] xi+n+1 −xi = = 4nfi+1 n!hn−4nfi n!hn (n+ 1)h=4nfi+1 − 4nfi (n+ 1)!hn+1 = =4n+1fi (n+ 1)!hn+1 , de donde se tiene el resultado.  Aplicando este lema e introduciendo la variable s= (x−x0)/h, podemos escribir (x−x0)(x−x1). . . (x−xk−1)f[x0, x1, . . . , xk] = s(s−1) . . . (s−k+ 1) k!4kf0,(1.25) Introducimos la notación: para s∈Rykentero,k≥1, s k=s(s−1) . . . (s−k+ 1) k!(1.26) Teniendo en cuenta (1.25) y (1.26), la fórmula de Newton para el polinomio de interpolación de frelativo a los n+1 puntos equidistantes x0, . . . , xnse puede simplificar como sigue Pn(x) = f0+s 14f0+s 242f0+· · · +s n4nf0(1.27) A esta forma de escribir el polinomio de interpolación se la suele denominar fórmula de Newton progresiva. De forma similar a lo que sucedía con las diferencias divididas, las diferencias progresivas pueden ser obtenidas mediante la Tabla 2. De este modo, una diferencia progresiva se obtiene restando dos diferencias de orden inferior. 11 Tabla 2. Diferencias finitas progresivas. xfi4fi42fi43fi. . . x0 x1 x2 x3 x4 . . . f0 f1 f2 f3 f4 . . . f1−f0=4f0 f2−f1=4f1 f3−f2=4f2 f4−f3=4f3 . . . 42f0 42f1 42f2 . . . 43f0 43f1 . . . . . . Lema 1.12. Suponemos que f(x)tiene derivada (n+ 1)-ésima en el intervalo [a, b]. Sea Pn(x)el polinomio de interpolación de f(x)con respecto a n+ 1 puntos equidistantes xi=x0+ih,i= 0,1, . . . , n pertenecientes al intervalo [a, b], siendo h > 0el paso. Entonces, para todo x=x0+sh perteneciente al intervalo [a, b]existe un punto ξ=ξ(x)en el intervalo abierto m´ın{x0, xn, x}< ξ < m´ax{x0, xn, x} tal que f(x)−Pn(x) = hn+1s n+ 1f(n+1)(ξ). Demostración: El resultado se obtiene directamente del Teorema 1.5. Bajo estas condiciones se tiene la representación del error (1.9); usando la notación (1.26) de esta sección y que s= (x−x0)/h se tiene el resultado.  Definición 1.13. Definimos el operador diferencia regresiva,O, por: O0fi=fi Ofi=fi−fi−1 12 y para todo entero k≥0 Ok+1fi=O(Okfi) = Okfi−Okfi−1. Lema 1.14. Se verifica que: f[xi−k, xi−k, . . . , xi] = Okfi k!hk,∀k≥0. Demostración: Análoga a la demostración del Lema 1.11.  Dado que Pn(x) = f[xn] + f[xn, xn−1](x−xn) + · · · +f[xn, . . . , x0](x−xn). . . (x−x1), haciendo el cambio de variable s= (x−xn)/h y usando la notación (1.26), obtenemos la expresión del polinomio de interpolación Pn(x) Pn(x) = fn+s 1Ofn+s+ 1 2O2fn+· · · +s+n−1 nOnfn,(1.28) que recibe el nombre de fórmula de Newton regresiva. Podemos calcularla fácilmente con el uso de la tabla de diferencias finitas regresivas (Tabla 3). 13 Ψn,j(xi) = 0,Ψ0 n,j(xi) = δij, i, j = 0,1, . . . , n. (2.3) Estos polinomios vienen dados, en términos de los polinomios fundamentales de Lagrange, φn,j(x), como: ψn,j(x)≡[1 −2φ0 n,j(xj)(x−xj)]φ2 n,j(x), Ψn,j(x)≡(x−xj)φ2 n,j(x).(2.4) Generalicemos este tipo de interpolación: Teorema 2.1. Dados k+ 1 puntos distintos x0, x1, . . . , xkde un intervalo [a, b]yk+ 1 enteros no negativos α0, α1, . . . , αk, denotamos n=k+α0+α1+· · · +αk. Dada una función fdefinida en el intervalo [a, b]que admite derivada de orden αien el punto xi para todo i= 0, . . . , k, existe un único polinomio Pnde grado menor o igual que n tal que: Para todo (i, l)tal que 0≤i≤k, 0≤l≤αi, se tiene que: P(l) n(xi) = f(l)(xi),(2.5) donde f(l)(xi)denota la derivada de orden lde fen el punto xi. Demostración: Las ecuaciones (2.5) constituyen un sistema lineal de n+1 ecuaciones y n+1 incógnitas (los coeficientes de Pn). Será suficiente demostrar que el sistema homogéneo asociado admite únicamente la solución nula, es decir que las relaciones Pn∈Pn,y∀(i, l)con 0≤i≤ky0≤l≤αi, P(l) n(xi)=0,(2.6) implican que Pn≡0. En efecto, (2.6) implica que para todo i= 0,1, . . . , k,xies una raíz de multiplicidad αi+ 1 de Pn; por tanto Pnes de la forma: Pn(x) = q(x) k Y i=0 (x−xi)αi+1, donde qes un polinomio. Como Pk i=0(αi+ 1) = n+ 1, esto sólo es compatible con el hecho de que Pn∈Pnsi q≡0, por tanto Pn≡0. Definición 2.2. El polinomio Pndefinido en el teorema anterior recibe el nombre de polinomio de interpolación de Hermite de la función frelativo a los puntos x0, x1, . . . , xk y a los enteros α0, α1, . . . , αk. Nótese que si tomamos αi= 0 para todo i= 0,1, . . . , k obtenemos la definición de polinomio de interpolación de Lagrange, mientras que si αi= 1 para todo i= 0,1, . . . , k llegamos al caso particular explicado al principio de este capítulo. 20 Cuando αi≥1para todo i= 0, . . . , k, se habla de interpolación osculatoria. En este caso, en cada punto xiel polinomio Pncoincide con fy la recta tangente a Pncoincide con la recta tangente a f. Además, si k= 0 yα0=nobtenemos el polinomio de Taylor de grado n. Supongamos dados los números reales bil para todo par (i, l)tal que 0≤i≤k, 0≤ l≤αi; en virtud del Teorema 2.1 el problema (Encontrar Pn∈Pntal que ∀(i, l)con 0≤i≤k, 0≤l≤αi, P(l) n(xi) = bil, tiene solución única. En particular, si para un par (i, l)tomamos bil = 1 ybjm = 0,∀(j, m)6= (i, l), obtenemos los polinomios de base o fundamentales,pil, para la interpolación de Hermite relativa a los puntos x0, x1, . . . , xky los enteros α0, α1, . . . , αk. Sea qi(x) = k Y j=0 j6=i x−xj xi−xjαj+1 , los polinomios de base pi,l están dados por: pi,αi(x) = (x−xi)αi αi!qi(x), y para l=αi−1, αi−2,...,1,0, pi,l(x) = (x−xi)l l!qi(x)− αi X j=l+1 j lq(j−l) i(xi)pi,j(x). Sea i= 0, . . . , k yl= 0, . . . , αi; tenemos que comprobar las condiciones: p(l) i,l (xi)=1, p(m) i,l (xi)=0,∀m= 0, . . . , αi, m 6=l p(m) i,l (xj) = 0,∀j= 0, . . . , k, j 6=iy∀m= 0, . . . , αj.(2.7) Procederemos de forma recursiva en l, en orden decreciente: l=αi, αi−1,...,0. Comenzamos pues con l=αi. Como pi,αi(x)contiene el factor (x−xj)αj+1 para todo j= 0, . . . , k,j6=i, tenemos que p(m) i,αi(xj) = 0,∀j= 0, . . . , k y∀m= 0, . . . , αj. 21 Dado que pi,αi(x)contiene el factor (x−xi)αi, tenemos que p(m) i,αi(xi)=0,∀m= 0, . . . , αi−1. Aplicando la regla de Leibniz, obtenemos p(αi) i,αi(x) = αi X t=0 αi tαi(αi−1) . . . (αi−y+ 1) αi!(x−xi)αi−tq(αi−t) i(x). Todos los sumandos se anulan en x=xi, salvo el correspondiente a t=αi, por tanto p(αi) i,αi(xi) = αi αiαi! αi!q(0) i(xi) = qi(xi) = 1. Consideramos ahora l=α1,...,0y suponemos ciertas las condiciones relativas a pi,r, para l+ 1 ≤r≤αi. pi,l(x)contiene el factor (x−xj)αj+1 para todo j= 0, . . . , k,j6=i, por tanto, p(m) i,l (xj) = 0,∀j= 0, . . . , k, j 6=i. Sea mentero, o≤m≤l−1, p(m) i,l (x) = dm dxm(x−xi)l l!qi(x)− αi X j=l+1 j lq(j−l) i(xi)p(m) i,j (x) El primer término del segundo miembro se anula en xi; por otra parte, para todo jcon αi≥j≥l+ 1, como j > m,p(m) i,j (xi)=0. Luego p(m) i,l (xi) = 0. Sea ahora m=l, p(l) i,l (x) = l X t=0 l tl(l−1) . . . (l−t+ 1) l!(x−xi)l−tq(l−t) i(x)− αi X j=l+1 j lq(j−l) i(xi)p(l) i,j(x). Todos los sumandos de la primera suma se anulan en xisalvo el correspondiente a t=l; por otra parte p(l) i,j(xi) = 0 ya que j > l. Obtenemos, por tanto p(l) i,l (xi) = l ll! l!q(0) i(xi)=1. Consideremos,por último, el caso en que l+ 1 ≤m≤αi, p(m) i,l (x) = l X t=0 m tl(l−1) . . . (l−t+ 1) l!(x−xi)l−tq(m−t) i(x)− αi X j=l+1 j lq(j−l) i(xi)p(m) i,j (x). 22 En la primera suma, el único término que no se anula en xies el que corresponde a t=ly su valor en dicho punto es m lq(m−l) i(xi). Por otra parte, los valores p(m) i,j (xi)son todos nulos salvo para j=m, ya que l+ 1 ≤ m≤αi. Por tanto, el valor de la segunda suma en xies m lq(m−l) i(xi)p(m) i,m (xi) = m lq(m−l) i(xi). Concluimos pues que p(m) i,l (xi)=0. Con esto quedan comprobadas todas las condiciones (2.7). El polinomio de interpolación de Hermite de fse obtiene a partir de los polinomios de base por: Pn(x) = k X i=0 αi X l=0 f(l)(xi)pil(x).(2.8) 2.2. Acotación del error para la interpolación de Hermite Teorema 2.3. Supongamos f∈Cn+1([a, b]). Para todo x∈[a, b], x 6=x0, . . . , xkexiste ξ=ξ(x)en el intervalo abierto: m´ın{x0, x1, ..., xk, x}< ξ < m´ax{x0, x1, ..., xk, x}(2.9) tal que E(x) = f(x)−Pn(x) = 1 (n+ 1)!πn(x)f(n+1)(ξ),(2.10) donde πn(x) = k Y i=0 (x−xi)αi+1. Demostración: Consideremos xfijo; definimos el polinomio de grado n+1 en la variable t qn+1(t) = Pn(t) + f(x)−Pn(x) πn(x)πn(t). Para todo (i, l)con 0≤i≤k, 0≤l≤αi,qn+1 verifica que q(l) n+1(xi) = f(l)(xi); además qn+1(x) = f(x). (qn+1 es el polinomio de interpolación de Hermite relativo a los puntos x0, x1, . . . , xk, x y a los enteros α0, α1, . . . , αk,0.) 23 Consideremos la función F(t) = f(t)−qn+1(t), que tiene en xiun cero de multiplicidad mayor o igual que αi+ 1 y en el punto xun cero de multiplicidad mayor o igual que 1. Así, el número de ceros de la función F, contados con sus multiplicidades, es al menos Pk i=0(αi+ 1) + 1 = n+ 2. En virtud del Teorema de Rolle, deducimos que F0tiene al menos n+1 ceros, contados con sus multiplicidades. En un número finito de pasos se llega a que la función F(n+1) admite al menos un cero, ξ, en el intervalo (2.9). Por tanto: 0 = F(n+1)(ξ) = f(n+1)(ξ)−q(n+1) n+1 (ξ).(2.11) El polinomio qn+1(t)es de grado n+ 1 y su coeficiente principal es f(x)−Pn(x) πn(x); en consecuencia, q(n+1) n+1 (t) = f(x)−Pn(x) πn(x)(n+ 1)!, lo que, junto con (2.11) implica (2.10).  Observación 2.4.Si x=xi,E(x)=0y (2.10) se verifica trivialmente. En conclusión, si podemos acotar |f(n+1)(x)|en [a, b], el teorema anterior nos proporciona una acotación del error. Ejemplo 2.5. Sea f(x) = tan(πx), supongamos que queremos calcular f(1/8) sabiendo que f(0) = 0,f(1/4) = 1,f0(0) = πyf0(1/4) = 2π. Tenemos que a=x0= 0,b=x1= 1/4yα0=α1= 1, con lo cual n= 3. Empleando las fórmulas (2.2) y (2.4) podemos calculamos el polinomio de interpolación de Hermite de grado menor o igual que 3. Empleando el programa HERMITE1.m, cuyo código se incluye en el anexo, se obtiene pol = 16*X^2 - 128*X^2*(X - 1/4) + 32*pi*X^2*(X - 1/4) + X*pi*(4*X - 1)^2. Haciendo x= 1/8obtenemos el valor aproximado de fen 1/8: P3(1/8) = 1 2−π 32 ≃0,4018 . . . Aplicando el Teorema 2.3 podemos acotar el error f(x)−P3(x) = f(4)(ξ) 4! x2(x−1/4)2 para algún ξ∈[0,1/4]. 24 f(4)(x)=8π4(2 tan xπ + 5 tan3xπ + 3 tan5xπ) En el intervalo [0,1/4],0≤tan xπ ≤1, por lo que, 0≤f(4)(x)≤80π4en dicho intervalo. Así, |f(1/8) −P3(1/8)| ≤ 5π4 3·211 ≃0,0792 . . . Definición 2.6. Sea f∈C([a, b]) definimos la norma ||.||C([a,b]) por ||f||C([a,b]) = sup x∈[a,b] |f(x)|. Veamos a continuación algunos ejemplos que ilustran el Teorema 2.3. Ejemplo 2.7. Sea f∈C4([a, b]), vamos a acotar el error para el polinomio de interpolación de Hermite P3, es decir, siendo k= 1,x0=a,x1=byα0=α1= 1, con lo que n= 3. P3interpola a f(x)y su primera derivada. Aplicamos el Teorema 2.3 para estimar f(x)−P3(x): f(x)−P3(x) = f(4)(ξ) 4! (x−a)2(x−b)2, con ξ∈(a, b). Como |(x−a)(x−b)| ≤ (b−a)2/4, se tiene |f(x)−P3(x)| ≤ ||f(4)||C([a,b]) 4!24(b−a)4. Ejemplo 2.8. Sea ahora f∈C6([a, b]), acotamos el error para el polinomio de interpolación de Hermite P5, siendo k= 1,x0=a,x1=byα0=α1= 2, con lo que n= 5. P5interpola a f(x)y las derivadas primera y segunda. Para estimar f(x)−P5(x)aplicamos el Teorema 2.3 y se tiene f(x)−P5(x) = f(6)(ξ) 6! (x−a)3(x−b)3, con ξ∈(a, b). Como |(x−a)(x−b)| ≤ (b−a)2/4se deduce que |f(x)−P5(x)| ≤ ||f(6)||C([a,b]) 6!26(b−a)6. Podemos obtener también una expresión para acotar f(r)−P(r) n, como se muestra en el siguiente teorema: 25 Teorema 2.9. Supongamos f∈Cn+1([a, b]). Sea Pnel polinomio de interpolación de Hermite de frespecto a los k+ 1 puntos x0, x1, . . . , xken el intervalo [a, b]y a los enteros α0, α1, . . . , αk, siendo n=k+α0+α1+· · · +αk. Sea 1≤r≤n.Entonces ||f(r)−P(r) n||C[a,b]≤||f(n+1)||C([a,b]) (n+ 1 −r)! (b−a)n+1−r(2.12) Demostración: Consideremos primero el caso k= 0, con lo que xk=x0yn=α0. El polinomio de interpolación en este caso es el polinomio de Taylor de grado nde fentorno al punto x0 Pn(x) = (Tnf)x0(x) = n X j=0 f(j)(x0) j!(x−x0)j Así, puesto que P(r) n(x)es el polinomio de Taylor de grado n−rde f(r)en torno a x0, por la representación del error del polinomio de Taylor se tiene |f(r)(x)−P(r) n(x)|=|f(r)(x)−(Tn−rf(r))x0(x)| ≤ ≤1 (n−r+ 1)! m´ax η∈I[x0,x]|f(n+1)(η)||x−x0|n−r+1,(2.13) con x, x0∈[a, b],I[a, b] = [m´ın(x0, x),m´ax(x0, x)]. Esta desigualdad implica (2.12). Consideremos ahora el caso k≥1.f(x)−P(x)tiene al menos n+ 1 ceros (contados con sus multiplicidades) en el intervalo J:= [m´ın(x0, . . . , xk),m´ax(x0, . . . , xk)]. Aplicando reiteradamente el Teorema de Rolle resulta que f(r)(x)−P(r) n(x)tiene al menos n+ 1 −r ceros (contados con sus multiplicidades) en el intervalo J. Nótese que n+ 1 −r≥1ya que r≤n. Por tanto, podemos encontrar en el intervalo J l ceros distintos de f(r)(x)−P(r) n(x), donde 1≤l≤n+ 1 −r, que denotamos ξ1, . . . , ξl, de manera que sus multiplicidades respectivas, ˆ β1,..., ˆ βl, cumplen ˆ β1+ˆ β2+· · · +ˆ βl≥n+ 1 −r. Sean β1, . . . , βlenteros tales que 1≤βj≤ˆ βjpara todo j= 1, . . . , l yβ1+· · · +βl=n+ 1 −r. Definimos la función auxiliar: E(t) = f(r)(t)−P(r) n(t)−f(r)(x)−P(r) n(x) Ql j=1(x−ξj)βj l Y j=1 (t−ξj)βj. Dado que f∈Cn+1([a, b]),E∈Cn+1−r([a, b]). Para cada j= 1, . . . , l,ξjes un cero de Econ multiplicidad mayor o igual que βj y como E(x) = 0,Etiene al menos β1+· · · +βl+ 1 = n−r+ 2 ceros, contados con sus multiplicidades, en el intervalo Jx:= [m´ın(x0, . . . , xk, x),m´ax(x0, . . . , xk, x)]. Por el Teorema de Rolle, E0tiene al menos n−r+ 1 ceros en el intervalo Jxcontados con sus multiplicidades. Aplicando el Teorema de Rolle repetidamente llegamos a que existe ξ∈Jx tal que E(n−r+1)(ξ)=0. 26 Dado que E(n−r+1)(t) = f(n+1)(t)−(n+ 1 −r)!f(r)(x)−P(r) n(x) Ql j=1(x−ξj)βj , se tiene f(r)(x)−P(r) n(x) = f(n+1)(ξ) (n+ 1 −r)! l Y j=1 (x−ξj)βj(2.14) Nótese que si x=ξjpara algún j= 1, . . . , l, esta ecuación también se cumple. Finalmente, (2.12) se deduce de (2.14).  Veamos algunos ejemplos. Ejemplo 2.10. Bajo las condiciones del Ejemplo 2.7, obtenemos las siguientes acotaciones: ||f0−P0||C([a,b]) ≤||f(4)||C([a,b]) 3! (b−a)3. ||f00 −P00||C([a,b]) ≤||f(4)||C([a,b]) 2(b−a)2. ||f(3) −P(3)||C([a,b]) ≤ ||f(4)||C([a,b])(b−a). Ejemplo 2.11. Bajo las condiciones del Ejemplo 2.8 podemos obtener las siguientes cotas: ||f0−P0||C([a,b]) ≤||f(6)||C([a,b]) 5! (b−a)5. ||f00 −P00||C([a,b]) ≤||f(6)||C([a,b]) 4! (b−a)4. ||f(3) −P(3)||C([a,b]) ≤||f(6)||C([a,b]) 3! (b−a)3. ||f(4) −P(4)||C([a,b]) ≤||f(6)||C([a,b]) 2(b−a)2. ||f(5) −P(5)||C([a,b]) ≤ ||f(6)||C([a,b])(b−a). Ejemplo 2.12. Sea f∈C2m([a, b]) yP2m−1el polinomio de interpolación, con k= 1, x0=a, x1=byα0=α1=m−1 P2m−1interpola a fy a sus derivadas hasta el orden m−1en ayb. Aplicando el Teorema 2.9, ||f(r)−P(r) 2m−1||C[a,b]≤||f(2m)||C([a,b]) (2m−r)! (b−a)2m−r,0≤r≤2m−1. Observación 2.13.Los ejemplos 2.10 y 2.11 son casos particulares del Ejemplo 2.12. En el capítulo 3de [6] se obtienen unas cotas similares a las dadas en los ejemplos 2.11 y 2.12. 27 28 Capítulo 3 Diferencias divididas con argumentos repetidos En este último capítulo nos centraremos en extender la definición de las diferencias divididas al caso de argumentos repetidos. Veremos también que esta nueva expresión de las diferencias divididas nos permitirá construir el polinomio de interpolación de Hermite mediante la elaboración de tablas como habíamos hecho con la fórmula de interpolación de Newton. Finalmente, estudiaremos el límite del polinomio de interpolación de Lagrange cuando los puntos de interpolación tienden a unos nodos límite, pudiendo coincidir algunos de estos. Seguiremos principalmente la sección 2.7 de [3]. 3.1. Diferencias divididas con argumentos repetidos Anteriormente definíamos f[x0, . . . , xk]solamente cuando todos los k+1 puntos x0, . . . , xk son distintos. Si tomamos k= 1 f[x0, x1] = f(x1)−f(x0) x1−x0 , que sólo tiene sentido si x16=x0. Si f(x)es continuamente diferenciable l´ım x1→x0 f(x1)−f(x0) x1−x0 =f0(x0), con lo que parece coherente definir la diferencia dividida f[x0, x1] = f0(x0) 29 la hipótesis de las repeticiones consecutivas. Esto es debido a que se usa la fórmula (3,8), obtenida en la demostración del Teorema 3.6, en la cual que usó dicha hipótesis. Veamos un ejemplo: Ejemplo 3.10. Vamos a construir la tabla de las diferencias divididas para calcular el polinomio de interpolación de f(x) = sinh xdados los valores: f(0) = 0, f0(0) = 1, f(1) = 1,1752, f0(1) = 1,5431 Obviamente, el polinomio de interpolación será de grado menor o igual que 3y usando (3.17) obtenemos la siguiente tabla: xf(x)f[x, x]f[x, x, x]f[x, x, x, x] 0 0 1 1 0 0 1,1752 1,1752 1 1,1752 1,5431 0,1752 0,3679 0,1927 ˆ P3(x) = 0 + 1(x−0) + (0,1752)(x−0)2+ (0,1927)(x−0)2(x−1) Si calculamos el valor aproximado de la función en 0,5se tiene que P3(0,5) = 0,5197125 mientras que f(0,5) = 0,5211. Puesto que los puntos de interpolación están bastante alejados entre si el error es considerable. El siguiente ejemplo está propuesto como ejercicio en [5]. Ejemplo 3.11. Consideramos los valores: sin(1,6) = 0,9995736030 cos(1,6) = −0,9291995223 sin(1,7) = 0,9916648105 cos(1,7) = −0,1288444943 36 Queremos aproximar el valor de sin(1,65). Tenemos que f(x) = sin(x),x0=a= 1,6y x1=b= 1,7, construyamos primero la tabla de las diferencias divididas para estos datos: xf(x)f[x, x]f[x, x, x]f[x, x, x, x] 1,6 1,6 1,7 1,7 0,9995736030 0,9995736030 0,9916648105 0,9916648105 −0,0291995223 −0,079087925 −0,1288444943 −0,498884027 −0,497565693 0,01318334 Se obtiene el polinomio de interpolación de Hermite: P_3(X) = (3799839047132805*X^3)/288230376151711744 - (324825484178482993*X^2)/576460752303423488 + (2410609107538067757*X)/1441151880758558720 - 259609249427398889/900719925474099200 Para los cálculos hemos empleado el programa HERMITE2.m, cuyo código se adjunta en el apéndice. Al evaluar este polinomio en el punto 1,65 se obtiene que P3(1,65) = 0,996864768900000 mientras que el valor real de la función es f(1,65) = 0,996865028453919. El error es por tanto |P3(1,65) −f(1,65)|= 2,59554 ·10−7(3.18) Acotemos el error cometido al aproximar f(x)por P3(x)en el intervalo [1,6,1,7], para ello usaremos el Teorema 2.9. f(4)(x) = sin(x)es decreciente y positiva en el intervalo [1,6,1,7], por tanto, ||f(4)||C([a,b]) =f(4)(1,6) = sin(1,6) = 0,9995736030, 37 y entonces |P3(x)−f(x)| ≤ 2,60306 ·10−7.(3.19) Vemos que el error calculado en (3.18) está por debajo de la cota hallada en (3.19). 38 ANEXO I. Codigo Matlab para el polinomio de interpolación de Lagrange Estas líneas de código nos permiten obtener la fórmula de Newton para el polinomio de interpolación de Lagrange. El programa NEWTON.m obtiene también una matriz triangular inferior que se corresponde con la tabla de las diferencias divididas. %POLINOMIO DE INTERPOLACIÓN DE NEWTON %ORDEN DEL POLINOMIO DE INTERPOLACIÓN n=input(’Introduca el grado del polinomio de interpolación’); %PUNTOS DE INTERPOLACIÓN x=zeros(1,n+1); fprintf(’Introduzca %d los puntos de interpolación: \n’,n+1) for i=1:n+1 x(i)=input(’’); end %VALORES DE LA FUNCIÓN fx=zeros(1,n+1); fprintf(’Introduzca los valores delafunción en dichos puntos’) for i=1:n+1 fx(i)=input(’’); end %CÁLCULO DE LA TABLA DE LAS DIFERENCIAS DIVIDIDAS. Obtendremos una matriz 39 %triangular inferior D=zeros(n+1,n+2); for i=1:n+1 D(i,1)=x(i); D(i,2)=fx(i); for j=3:n+2 if i+1>=j D(i,j)=(D(i,j-1)-D(i-1,j-1))/(D(i,1)-D(i-1-j+3,1)); end end end fprintf(’La matriz de las diferencias divididas es:\n’) disp(D) %CONSTRUIMOS EL POLINOMIO DE INTERPOLACIÓN syms X Q=D(1,2); p=0; for i=2:n+1 p=D(i,i+1); for j=2:i; p=p*(X-D(j-1,1)); end Q=Q+p; end fprintf(’El polinomio de interpolación de Newton es’) disp(Q) %GRÁFICA %Representamos ahora el polinomio de interpolación en azul. Los puntos de %interpolación están señalados con asteriscos rojos. a=input(’Introduce el exremo inferior de la gráfica’); b=input(’Introduce el extremo superior de la gráfica’); T=linspace(a,b); Y=subs(Q,X,T); 40 plot(T,Y) hold on plot(x,fx,’r*’) %CÁLCULO APROXIMADO DEL VALOR DE LA FUNCIÓN EN UN PUNTO fprintf(’Introduce un punto donde quieras aproximar el valor de la función’) t=input(’’); subs(Q,t) 41 42 ANEXO II. Codigo Matlab para el polinomio de interpolación de Hermite Hemos visto a lo largo de este trabajo que podemos obtener el polinomio de Hermite de dos maneras distintas: mediante la fórmula (2.8) y mediante diferencias divididas. Se muestran a continuación las líneas de código del programa HERMITE1.m que obtiene el polinomio de interpolación de Hermite usando la fórmula (2.8). También se incluyen los datos correspondientes al Ejemplo 2.5 en el archivo data1.m. %CÁLCULO DEL POLINOMIO DE HERMITE USANDO LA FÓRMULA (2.8) clc clear all data2 k=length(x)-1; m=[fx’,dx]; syms X for i=1:k+1 prod=1; for j=1:k+1 if i~=j res=(((X-x(j))/(x(i)-x(j)))^(alfa(j)+1)); prod=prod*res; end 43 end q(i)=prod; end for i=1:k+1 syms X p(i,alfa(i)+1)=((X-x(i))^alfa(i)/factorial(alfa(i)))*q(i); for l=alfa(i):-1:1 p(i,l)=((X-x(i))^(l-1)/factorial(l-1))*q(i); sum=0; for j=l:alfa(i) dq(i)=diff(q(i),j-l+1); qixi=subs(dq(i),x(i)); sum=sum+nchoosek(j,l-1)*qixi*p(i,j+1); end p(i,l)=p(i,l)-sum; end end pol=0; for i=1:k+1 for l=1:alfa(i)+1 pol=m(i,l)*p(i,l)+pol; end end pol Código de data1.m relativo al Ejemplo 2.5: 44 %DATOS DEL EJEMPLO 2.5 %n= grado del polinomio de interpolación %x= puntos de interpolación %fx valores de la funcion %dx matriz de derivadas x=[0 0.25] fx=[0 1] dx=[pi;2*pi] alfa=[1 1] A continuación se muestran las líneas de código del programa HERMITE2.m para obtener el polinomio de interpolación de Hermite en la forma de Newton, así como la tabla de las diferencias divididas en forma de matriz triangular inferior. Se dan los datos del Ejemplo 3.11 en el archivo data2.m, mostrado abajo. %CÁLCULO DEL POLINOMIO DE HERMITE EN LA FORMA DE NEWTON clc clear all data D=[x’,fx’,dx] DD=zeros(n+1,n+2); %MATRIZ DE LAS DIFERENCIAS DIVIDIDAS for i=1:n+1 DD(i,1)=x(i); DD(i,2)=fx(i); end for i=2:n+1 for j=3:n+2 if i+1>=j if D(i,1)==D(i-1-j+3,1) DD(i,j)=D(i-1,j)/factorial(j-2); else DD(i,j)=(DD(i,j-1)-DD(i-1,j-1))/(DD(i,1)-DD(i-1-j+3,1)); 45