Burden Issuu.pdf 6c713j

10a. edición

ANÁLISIS

NUMÉRICO Richard L. Burden • Douglas J. Faires • Annette M. Burden

Análisis numérico DÉCIMA EDICIÓN

Richard L. Burden Youngstown University

J. Douglas Faires Youngstown University

Annette M. Burden Youngstown University

Traducción: Mara Paulina Suárez Moreno Traductora profesional

Revisión técnica: Wilmar Alberto Díaz Ossa Mágister en matemáticas aplicadas Profesor en la Universidad Distrital Francisco José de Caldas

Australia • Brasil • Corea • España • Estados Unidos • Japón • México • Reino Unido • Singapur

Análisis numérico, 10 a. ed. Richard L. Burden, J. Douglas Faires y Annette M. Burden Director Editorial para Latinoamérica: Ricardo H. Rodríguez Editora de Adquisiciones para Latinoamérica: Claudia C. Garay Castro Gerente de Manufactura para Latinoamérica: Antonio Mateos Martínez Gerente Editorial de Contenidos en Español: Pilar Hernández Santamarina Gerente de Proyectos Especiales: Luciana Rabuffetti Coordinador de Manufactura: Rafael Pérez González Editora: Ivonne Arciniega Torres Diseño de portada: Anneli Daniela Torres Arroyo Imagen de portada: © theromb/Shutterstock.com Composición tipográfica: Tsuki Marketing S.A. de C.V. Gerardo Larios García

Impreso en México 1 2 3 4 5 6 7 20 19 18 17

© D.R. 2017 por Cengage Learning Editores, S.A. de C.V., una Compañía de Cengage Learning, Inc. Corporativo Santa Fe Av. Santa Fe núm. 505, piso 12 Col. Cruz Manca, Santa Fe C.P. 05349, México, D.F. Cengage Learning® es una marca registrada usada bajo permiso. DERECHOS RESERVADOS. Ninguna parte de este trabajo amparado por la Ley Federal del Derecho de Autor, podrá ser reproducida, transmitida, almacenada o utilizada en cualquier forma o por cualquier medio, ya sea gráfico, electrónico o mecánico, incluyendo, pero sin limitarse a lo siguiente: fotocopiado, reproducción, escaneo, digitalización, grabación en audio, distribución en internet, distribución en redes de información o almacenamiento y recopilación en sistemas de información a excepción de lo permitido en el Capítulo III, Artículo 27 de la Ley Federal del Derecho de Autor, sin el consentimiento por escrito de la Editorial. Traducido del libro Numerical Analysis, Tenth Edition Richard L. Burden, J. Douglas Faires, Annette M. Burden Publicado en inglés por Cengage Learning © 2016, 2011, 2005 ISBN: 978-1-305-25366-7 Datos para catalogación bibliográfica: Burden, Faires y Burden Análisis numérico, 10a. ed. ISBN: 978-607-526-411-0 Visite nuestro sitio en: http://latinoamerica.cengage.com

Contenido Prefacio

1

2

Preliminares matemáticos y análisis de error 1 1.1 1.2

Revisión de cálculo 2 Errores de redondeo y aritmética computacional 11

1.3

Algoritmos y convergencia

1.4

Software numérico

28

El método de bisección 36 ,WHUDFLyQGHSXQWRÀMR Método de Newton y sus extensiones 49 Análisis de error para métodos iterativos 58 Convergencia acelerada 64 Ceros de polinomios y método de Müller 68 Software numérico y revisión del capítulo 76

Interpolación y aproximación polinomial 3.1 3.2 3.3 3.4 3.5 3.6 3.7

4

22

Soluciones de las ecuaciones en una variable 35 2.1 2.3 2.4 2.5 2.6 2.7

3

vii

Interpolación y el polinomio de Lagrange 78 Aproximación de datos y método de Neville 86 Diferencias divididas 91 Interpolación de Hermite 99 Interpolación de spline cúbico 105 Curvas paramétricas 121 Software numérico y revisión del capítulo 126

Diferenciación numérica e integración 4.1 4.2 4.3

77

127

Diferenciación numérica 128 Extrapolación de Richardson 136 Elementos de integración numérica 142 iii

iv

Contenido

4.4 4.5 4.6 4.7 4.8 4.9 4.10

5

Problemas de valor inicial para ecuaciones de diferenciales ordinarias 193 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12

6

Teoría elemental de problemas de valor inicial 194 Método de Euler 198 Métodos de Taylor de orden superior 205 Método Runge-Kutta 209 Control de error y método Runge-Kutta-Fehlberg 218 Métodos multipasos 224 Método multipasos de tamaño de paso variable 236 Métodos de extrapolación 241 Ecuaciones de orden superior y sistemas de ecuaciones diferenciales 247 Estabilidad 254 Ecuaciones diferenciales rígidas 262 Software numérico 268

Métodos directos para resolver sistemas lineales 6.1 6.2 6.3 6.4 6.5 6.6 6.7

7

Integración numérica compuesta 150 Integración de Romberg 156 Métodos de cuadratura adaptable 162 Cuadratura gaussiana 168 Integrales múltiples 174 Integrales impropias 186 Software numérico y revisión del capítulo 191

Sistemas de ecuaciones lineales 270 Estrategias de pivoteo 279 Álgebra lineal e inversión de matriz 287 Determinante de una matriz 296 Factorización de matriz 298 Tipos especiales de matrices 306 Software numérico 318

Técnicas iterativas en álgebra de matrices 319 7.1 7.2 7.3 7.7

Normas de vectores y matrices 320 Eigenvalores y eigenvectores 329 Técnicas iterativas de Jacobi y Gauss-Siedel 334 7pFQLFDVGHUHODMDFLyQSDUDUHVROYHUVLVWHPDVOLQHDOHV &RWDVGHHUURU\UHÀQDPLHQWRLWHUDWLYR (OPpWRGRGHJUDGLHQWHFRQMXJDGR Software numérico 366

269

Contenido

8

Teoría de aproximación 369 8.1 8.2 8.3 8.4 8.5 8.6 8.7

9

Aproximación de eigenvalores 421 9.1 9.2 9.3 9.4 9.5 9.6 9.7

10

Aproximación por mínimos cuadrados discretos 370 Polinomios ortogonales y aproximación por mínimos cuadrados 378 Polinomios de Chebyshev y ahorro de series de potencia 385 Aproximación de función racional 393 Aproximación polinomial trigonométrica 402 Transformadas rápidas de Fourier 410 Software numérico 419

Álgebra lineal y eigenvalores 422 Matrices ortogonales y transformaciones de similitud 428 El método de potencia 431 Método de Householder 445 El algoritmo QR 452 Descomposición en valores singulares 462 Software numérico 474

Soluciones numéricas de sistemas de ecuaciones no lineales 475 3XQWRVÀMRVSDUDIXQFLRQHVGHYDULDVYDULDEOHV 10.2 Método de Newton 482 10.3 Métodos cuasi-Newton 487 10.4 Técnicas de descenso más rápido 492 10.5 Homotopía y métodos de continuación 498 10.6 Software numérico 504

11

Problemas de valor en la frontera para ecuaciones diferenciales ordinarias 505 11.1 El método de disparo lineal 506 11.2 El método de disparo para problemas no lineales 512 0pWRGRVGHGLIHUHQFLDVÀQLWDVSDUDSUREOHPDVOLQHDOHV 0pWRGRVGHGLIHUHQFLDVÀQLWDVSDUDSUREOHPDVOLQHDOHV 11.5 El método de Rayleigh-Ritz 527 11.6 Software numérico 540

v

vi

Contenido

12

Soluciones numéricas para ecuaciones diferenciales parciales 541 12.1 Ecuaciones diferenciales parciales elípticas 544 12.2 Ecuaciones diferenciales parciales parabólicas 551 12.3 Ecuaciones diferenciales parciales hiperbólicas 562 8QDLQWURGXFFLyQDOPpWRGRGHHOHPHQWRVÀQLWRV 12.5 Software numérico 579

Material en línea El siguiente material se encuentra disponible en línea: • • • • • • • • • • •

&RQMXQWRVGHHMHUFLFLRV Preguntas de análisis Conceptos clave Revisión de capítulo Bibliografía 5HVSXHVWDVDHMHUFLFLRVVHOHFFLRQDGRV Índice Índice de algoritmos Glosario de notación Trigonometría *UiÀFDVFRPXQHV

Ingrese a www.cengage.com, busque el libro por el ISBN e ingrese el siguiente código de :

CAPÍTULO

1

Preliminares matemáticos y análisis de error Introducción Al comenzar los cursos de química, estudiamos la ley del gas ideal, PV = NRT, que relaciona la presión P, el volumen V, la temperatura T y el número de moles N de un gas “ideal”. En esta ecuación, R es una contante que depende del sistema de medición. Suponga que se realizan dos experimentos para evaluar esta ley, mediante el mismo gas en cada caso. En el primer experimento,

P = 1.00 atm,

V = 0.100 m3 ,

N = 0.00420 mol,

R = 0.08206.

La ley del gas ideal predice que la temperatura del gas es

T =

(1.00)(0.100) PV = = 290.15 K = 17◦ C. NR (0.00420)(0.08206)

Sin embargo, cuando medimos la temperatura del gas, encontramos que la verdadera temperatura es 15◦C.

V1 V2

A continuación, repetimos el experimento utilizando los mismos valores de R y N, pero incrementamos la presión en un factor de dos y reducimos el volumen en ese mismo factor. El producto PV sigue siendo el mismo, por lo que la temperatura prevista sigue siendo 17◦C. Sin embargo, ahora encontramos que la temperatura real del gas es 19◦C. 1

2

CAPÍTULO 1

Preliminares matemáticos y análisis de error

Claramente, se sospecha la ley de gas ideal, pero antes de concluir que la ley es inválida en esta situación, deberíamos examinar los datos para observar si el error se puede atribuir a los resultados del experimento. En este caso, podríamos determinar qué tan precisos deberían ser nuestros resultados experimentales para evitar que se presente un error de esta magnitud. El análisis del error involucrado en los cálculos es un tema importante en análisis numérico y se presenta en la sección 1.2. Esta aplicación particular se considera en el ejercicio 26 de esa sección. Este capítulo contiene una revisión breve de los temas del cálculo de una sola variable que se necesitarán en capítulos posteriores. Un conocimiento sólido de cálculo es fundamental para comprender el análisis de las técnicas numéricas y sería preciso efectuar una revisión más rigurosa para quienes no han estado en o con este tema durante un tiempo. Además, existe una introducción a la convergencia, el análisis de error, la representación GHQ~PHURVHQOHQJXDMHGHPiTXLQD\DOJXQDVWpFQLFDVSDUDFODVLÀFDU\PLQLPL]DUHOHUURU computacional.

1.1 Revisión de cálculo Límites y continuidad Los conceptos de límite y continuidad de una función son fundamentales para el estudio del cálculo y constituyen la base para el análisis de las técnicas numéricas. Definición 1.1

Una función fGHÀQLGDHQXQFRQMXQWRX de números reales que tiene el límite L a x0, escrita como lím f (x) = L , x→x0

si, dado cualquier número real ε > 0, existe un número real δ > 0, de tal forma que

| f (x) − L| < ε,

siempre que

x∈X y

0 < |x − x0 | < δ.

FRQVXOWHODÀJXUD

Figura 1.1

ε

y

y 5 f (x)

L 1e L L 2e

x0 2 d

x0

x0 1 d

x

1.1

Definición 1.2 Los conceptos básicos de cálculo y sus aplicaciones se GHVDUUROODURQDÀQDOHVGHOVLJOR XVII y a principios del XVIII, pero los conceptos matemáticamente precisos de límites y continuidad se describieron hasta la época de Augustin Louis Cauchy ² +HLQULFK(GXDUG +HLQH² \.DUO :HLHUVWUDVV² DÀQDOHV del siglo XIX.

Definición 1.3

Revisión de cálculo

3

Sea f XQDIXQFLyQGHÀQLGDHQXQFRQMXQWRX de números reales y x0 ∈ X. Entonces f es continua en x0 si lím f (x) = f (x0 ). x→x0

La función f es continua en el conjunto X si es continua en cada número en X. El conjunto de todas las funciones que son continuas en el conjunto X se denota como C(X &XDQGRX es un intervalo de la recta real, se omiten los paréntesis en esta notación. Por ejemplo, el conjunto de todas las funciones continuas en el intervalo cerrado [a, b] se denota como C[a, b]. El símbolo R denota el conjunto de todos los números reales, que también tiene la notación del intervalo (−∞, ∞ 3RUHVRHOFRQMXQWRGHWRGDVODVIXQFLRQHVTXHVRQ continuas en cada número real se denota mediante C(R o mediante C(−∞, ∞ El límite de una sucesiónGHQ~PHURVUHDOHVRFRPSOHMRVVHGHÀQHGHPDQHUDVLPLODU Sea {xn }∞ n=1XQDVXFHVLyQLQÀQLWDGHQ~PHURVUHDOHV(VWDVXFHVLyQWLHQHHOlímite x (converge a x VLSDUDFXDOTXLHU ε 0, existe un entero positivo N(ε) tal que |xn − x| < ε siempre que n > N(ε /DQRWDFLyQ

lím xn = x,

n→∞

o

xn → x

en

n → ∞,

VLJQLÀFDTXHODVXFHVLyQ{xn }∞ n=1 converge a x. Teorema 1.4

Si fHVXQDIXQFLyQGHÀQLGDHQXQFRQMXQWRX de números reales y x0 ∈ X, entonces los siguientes enunciados son equivalentes: a.

f es continua en x0;

b.

Si {xn }∞ n=1 es cualquier sucesión en X, que converge a x0, entonces lím n→∞ f (xn ) = f (x0 ).

Se asumirá que las funciones que consideraremos al analizar los métodos numéricos son continuas porque éste es el requisito mínimo para una conducta predecible. Las funciones TXHQRVRQFRQWLQXDVSXHGHQSDVDUSRUDOWRSXQWRVGHLQWHUpVORFXDOSXHGHFDXVDUGLÀFXOWDdes al intentar aproximar la solución de un problema.

Diferenciabilidad /DVVXSRVLFLRQHVPiVVRÀVWLFDGDVVREUHXQDIXQFLyQSRUORJHQHUDOFRQGXFHQDPHMRUHVUHVXOWDGRVGHDSUR[LPDFLyQ3RUHMHPSORQRUPDOPHQWHXQDIXQFLyQFRQXQDJUiÀFDVXDYHVH comportaría de forma más predecible que una con numerosas características irregulares. La condición de uniformidad depende del concepto de la derivada. Definición 1.5

Si f HVXQDIXQFLyQGHÀQLGDHQXQLQWHUYDORDELHUWRTXHFRQWLHQHx0. La función f es diferenciable en x0 si f (x) − f (x0 ) f (x0 ) = lím x→x0 x − x0 existe. El número f (x0 ) recibe el nombre de derivada de f en x0. Una función que tiene una derivada en cada número en un conjunto X es diferenciable en X. La derivada de f en x0HVODSHQGLHQWHGHODUHFWDWDQJHQWHDODJUiÀFDGHf en (x0 , f (x0 )), FRPRVHPXHVWUDHQODÀJXUD

4

CAPÍTULO 1

Preliminares matemáticos y análisis de error

Figura 1.2 y La recta tangente tiene una pendiente f 9(x0)

f (x 0)

(x 0, f (x 0))

y 5 f (x)

x0

Teorema 1.6 El teorema atribuido a Michel 5ROOH² DSDUHFLy en 1691 en un tratado poco conocido titulado Méthode pour résoundre les égalites (Método para resolver las igualdades). Originalmente, Rolle criticaba el cálculo desarrollado por Isaac Newton y Gottfried Leibniz, pero después se convirtió en uno de sus defensores.

Teorema 1.7

x

Si la función f es diferenciable en x0, entonces f es continua en x0. Los siguientes teoremas son de importancia fundamental al deducir los métodos para estimación del cálculo de error. Las pruebas de estos teoremas y los otros resultados sin referencias en esta sección se pueden encontrar en cualquier texto de cálculo estándar. El conjunto de todas las funciones que tienen derivadas continuas n en X se denota como Cn(X y el conjunto de funciones que tienen derivadas de todos los órdenes en X se denota como C ∞ (X ). Las funciones polinomial, racional, trigonométrica, exponencial y logarítmica se encuentran en C ∞ (X ), donde X FRQVLVWHHQWRGRVORVQ~PHURVSDUDORVTXHVHGHÀQHQ las funciones. Cuando X es un intervalo de la recta real, de nuevo se omiten los paréntesis en esta notación. (Teorema de Rolle) Suponga que f ∈ C[a, b] y f es diferenciable en (a, b 6Lf(a = f (b HQWRQFHVH[LVWHXQQ~mero c en (a, b) con f 9(c = &RQVXOWHODÀJXUD

Figura 1.3

y

f 9(c) 5 0 y 5 f (x)

f (a) 5 f(b)

a

Teorema 1.8

c

b

x

(Teorema del valor medio) Si f ∈ C[a, b] y f es diferenciable en (a, b HQWRQFHV H[LVWH XQ Q~PHUR c en (a, b con FRQVXOWHODÀJXUD

f (c) =

f (b) − f (a) . b−a

1.1

Revisión de cálculo

5

Figura 1.4 y Líneas paralelas Pendiente f 9(c)

Pendiente

f (b) 2 f (a) b2a

c

a

Teorema 1.9

y 5 f (x)

x

b

(Teorema del valor extremo) Si f ∈ C[a, b], entonces existe c1, c2 ∈ [a, b] con f (c1 ≤ f (x ≤ f (c2 SDUDWRGDVODVx ∈ [a, b]. Además, si f es diferenciable en (a, b HQWRQFHVVHSUHVHQWDQORVQ~PHURVc1 y c2 ya sea en los extremos de [a, b] o donde f 9HVFHUR&RQVXOWHODÀJXUD

Figura 1.5 y

y 5 f (x)

c2

a

Ejemplo 1

c1

b

x

Encuentre los valores mínimo absoluto y máximo absoluto de f (x = 2 − ex + 2x en los intervalos a >, 1] y b >, 2]. Solución

Comenzamos por derivar f (x SDUDREWHQHU f v(x = −ex + 2.

f v(x = 0 cuando −ex + 2 = 0 o de forma equivalente, cuando ex = 2. Al tomar el logaritmo natural de ambos lados de la ecuación obtenemos ln (ex =OQ Rx =OQ ≈

6

CAPÍTULO 1

Preliminares matemáticos y análisis de error

a)

Cuando el intervalo es [0, 1], el extremo absoluto debe ocurrir en f , f OQ R f $OHYDOXDUWHQHPRV

f (0) = 2 − e0 + 2(0) = 1 f (ln (2)) = 2 − eln (2) + 2 ln (2) = 2 ln (2) ≈ 1.38629436112 f (1) = 2 − e + 2(1) = 4 − e ≈ 1.28171817154. Por lo tanto, el mínimo absoluto de f (x HQ > @ HV f = 1 y el máximo absoluto es f OQ =OQ b)

Cuando el intervalo es [1, 2], sabemos que f (x) = 0, por lo que el extremo absoluto se presenta en f \f 3RUORWDQWR f (2) = 2 − e2 + 2(2) = 6 − e2 ≈ . −1.3890560983 El mínimo absoluto en [1, 2] es 6 − e2 y el máximo absoluto es 1. Observamos que

máx | f (x)| = |6 − e2 | ≈ 1.3890560983.

0≤x≤2

En general, el siguiente teorema no se presenta en un curso de cálculo básico, pero se deriva al aplicar el teorema de Rolle sucesivamente a f, f , . . . ,\ÀQDOPHQWHD f (n−1) . Este resultado se considera en el ejercicio 26. Teorema 1.10

(Teorema generalizado de Rolle) Suponga que f ∈ C[a, b] es n veces diferenciable en (a, b 6Lf (x = 0 en los n + 1 números distintos a a ≤ x0 < x1 < 7 < xn ≤ b, entonces un número c en (x0, xn \SRUORWDQWRHQa, b existe con f (n (c = 0. También utilizaremos con frecuencia el teorema del valor intermedio. A pesar de que esta declaración parece razonable, su prueba va más allá del alcance del curso habitual de cálculo. Sin embargo, se puede encontrar en muchos textos de análisis (consulte, por ejemSOR>)X@S

Teorema 1.11

(Teorema del valor intermedio) Si f ∈ C[a, b] y K es cualquier número entre f(a \f (b HQWRQFHVH[LVWHXQQ~PHURc en (a, b SDUDHOFXDOf (c = K. /DÀJXUDPXHVWUDXQDRSFLyQSDUDHOQ~PHURJDUDQWL]DGDSRUHOWHRUHPDGHOYDORU intermedio. En este ejemplo, existen otras dos posibilidades.

Figura 1.6 y f (a)

(a, f (a)) y 5 f (x)

K f (b)

(b, f (b)) a

c

b

x

1.1

Ejemplo 2

Revisión de cálculo

7

Muestre que x5 − 2x +x2 − 1 = 0 tiene una solución en el intervalo [0, 1]. Solución &RQVLGHUH OD IXQFLyQ GHÀQLGD SRU f (x = x5 − 2x + x2 − 1. La función f es

continua en [0, 1]. Además, f = −1 < 0

y

0 < 1 = f .

Por lo tanto, el teorema del valor intermedio implica que existe un número c, con 0 , c , 1, para el cual c5 − 2c +c2 − 1 = 0. Como se observa en el ejemplo 2, el teorema del valor intermedio se utiliza para determinar cuándo existen soluciones para ciertos problemas. Sin embargo, no provee un medio HÀFLHQWHSDUDHQFRQWUDUHVWDVVROXFLRQHV(VWHWHPDVHFRQVLGHUDHQHOFDStWXOR

Integración El otro concepto básico del cálculo que se utilizará ampliamente es la integral de Riemann. Definición 1.12 George Fredrich Berhard Riemann ² UHDOL]yPXFKRVGH los descubrimientos importantes SDUDFODVLÀFDUODVIXQFLRQHV que tienen integrales. También realizó trabajos fundamentales en geometría y la teoría de funciones complejas y se le considera uno de los matemáticos prolíferos del siglo XIX.

La integral de Riemann de la función f en el intervalo [a, b] es el siguiente límite, siempre y cuando exista: b a

n

f (x) d x =

lím

máx xi →0

f (z i

xi ,

i=1

donde los números x0, x1,7 , xn satisfacen a = x0 ≤ x1 ≤ 7 ≤ xn = b, donde 6xi = xi − xi−1, para cada i = 1, 2,7 , n, y zi se selecciona de manera arbitraria en el intervalo [ xi−1 , xi ]. Una función f que es continua en un intervalo [a, b] es también Riemann integrable en [a, b]. Esto nos permite elegir, para conveniencia computacional, los puntos xi se separarán equitativamente en [a, b] para cada i = 1, 2, 7 , n, para seleccionar zi = xi. En este caso, b a

b−a f (x) d x = lím n→∞ n

n

f (xi ), i=1

GRQGHORVQ~PHURVPRVWUDGRVHQODÀJXUDFRPRxi, son xi = a + i(b − a n.

Figura 1.7 y y 5 f (x)

a 5 x 0 x1

x2 . . . x i21 x i

...

x n21 b 5 x n

x

Se necesitarán otros dos resultados en nuestro estudio para análisis numérico. El primero es una generalización del teorema del valor promedio para integrales.

8

CAPÍTULO 1

Preliminares matemáticos y análisis de error

Teorema 1.13

(Teorema del valor promedio para integrales) Suponga que f ∈ C[a, b], la integral de Riemann de g existe en [a, b], y g(x QRFDPELDGH signo en [a, b]. Entonces existe un número c en (a, b FRQ b a

f (x)g(x) d x = f (c)

b a

g(x) d x.

Cuando g(x ≡HOWHRUHPDHVHOWHRUHPDGHOYDORUPHGLRSDUDLQWHJUDOHVeVWH proporciona el valor promedio de la función f sobre el intervalo [a, b] como (consulte la ÀJXUD

f (c) =

1 b−a

b a

f (x) d x.

Figura 1.8 y y 5 f (x) f (c)

a

c

b

x

(QJHQHUDOODSUXHEDGHOWHRUHPDQRVHGDHQXQFXUVREiVLFRGHFiOFXORSHURVH SXHGHHQFRQWUDUHQPXFKRVWH[WRVGHDQiOLVLVFRQVXOWHSRUHMHPSOR>)X@S

Polinomios y series de Taylor (OWHRUHPDÀQDOHQHVWDUHYLVLyQGHFiOFXORGHVFULEHORVSROLQRPLRVGH7D\ORU(VWRVSROLQRmios se usan ampliamente en el análisis numérico. Teorema 1.14 %URRN7D\ORU² describió esta serie en 1715 en el artículo Methodus incrementorum directa et inversa (Métodos para incrementos directos e inversos). Isaac Newton, James Gregory y otros ya conocían algunos casos especiales del resultado y, probablemente, el resultado mismo.

(Teorema de Taylor) Suponga que f ∈ C n [a, b], f (n + existe en [a, b], y x0 ∈ [a, b]. Para cada x ∈ [a, b], existe un número ξ(x) entre x0 y x con

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

Pn (x) = f (x0 ) + f (x0 )(x − x0 ) + n

= k=0

f (k) (x0 ) (x − x0 )k k!

f (x0 ) f (n) (x0 ) (x − x0 )2 + · · · + (x − x0 )n 2! n!

1.1

Revisión de cálculo

9

y

Rn (x) = &ROLQ0DFODXULQ HV más conocido como el defensor del cálculo de Newton cuando éste fue objeto de los ataques LPSODFDEOHVGHORELVSR\ÀOyVRIR irlandés George Berkeley. Maclaurin no descubrió la serie que lleva su nombre; los matemáticos del siglo ya la conocían desde antes de que él naciera. Sin embargo, concibió un método para resolver un sistema de ecuaciones lineales que se conoce como regla de Cramer, que Cramer no publicó hasta 1750.

Ejemplo 3

f (n+1) (ξ(x)) (x − x0 )n+1 . (n + 1)!

Aquí Pn(x HV OODPDGR HO n-ésimo polinomio de Taylor para f alrededor de x0 y Rn(x recibe el nombre de residuo (o error de truncamiento UHODFLRQDGRFRQPn(x 3XHVWRTXHHO número ξ(x) en el error de truncamiento Rn(x GHSHQGHGHOYDORUGHx donde se evalúa el polinomio Pn(x HVXQDIXQFLyQGHODYDULDEOHx. Sin embargo, no deberíamos esperar ser capaces de determinar la función ξ(x) de manera explícita. El teorema de Taylor simplemente garantiza que esta función existe y que su valor se encuentra entre x y x0. De hecho, uno de los problemas comunes en los métodos numéricos es tratar de determinar un límite realista para el valor de f (n+1) (ξ(x)) cuando xVHHQFXHQWUDHQXQLQWHUYDORHVSHFtÀFR /DVHULHLQÀQLWDREWHQLGDDOWRPDUHOOtPLWHGHPn(x FRQIRUPH n → ∞ recibe el nombre de serie de Taylor para f alrededor de x0. En caso de que x0 = 0, entonces al polinomio de Taylor con frecuencia se le llama polinomio de Maclaurin y a la serie de Taylor a menudo se le conoce como serie de Maclaurin. El término error de truncamientoHQHOSROLQRPLRGH7D\ORUVHUHÀHUHDOHUURULPSOLFDGRDOXWLOL]DUXQDVXPDWUXQFDGDRÀQLWDSDUDDSUR[LPDUODVXPDGHXQDVHULHLQÀQLWD Si f(x = cos x y x0 = 0. Determine a)

el segundo polinomio de Taylor para f alrededor de x0; y

b)

el tercer polinomio de Taylor para f alrededor de x0.

Solución

Puesto que f ∈ C ∞ (R), el teorema de Taylor puede aplicarse a cualquiera n ≥ 0.

Además,

f (x) = − sen x, f (x) = − cos x, f (x) = sen x, y

f (4) (x) = cos x,

Por lo tanto f (0) = 1, f (0) = 0, f (0) = −1, a)

y

f (0) = 0.

Para n = 2 y x0 = 0, obtenemos cos x = f (0) + f (0)x +

f (0) 2 f (ξ(x)) 3 x + x 2! 3!

1 1 = 1 − x 2 + x 3 sen ξ(x), 2 6 donde ξ(x)HVDOJ~QQ~PHURSRUORJHQHUDOGHVFRQRFLGR HQWUH\x&RQVXOWHODÀJXUD Figura 1.9 y

1 p 22 2

2p

y 5 cos x p 2 2 p

1 2

y 5 P2(x) 5 1 2 2 x 2

x

10

CAPÍTULO 1

Preliminares matemáticos y análisis de error

Cuando x = 0.01, esto se convierte en

1 1 10−6 cos 0.01 = 1 − (0.01)2 + (0.01)3 sen ξ(0.01) = 0.99995 + sen ξ(0.01). 2 6 6 Por lo tanto, la aproximación para cos 0.01 provista por el polinomio de Taylor es 0.99995. El error de truncamiento, o término restante, relacionado con esta aproximación es

10−6 sen ξ(0.01) = 0.16 × 10−6 sen ξ(0.01), 6 donde la barra sobre el 6 en 0 .16 VHXWLOL]DSDUDLQGLFDUTXHHVWHGtJLWRVHUHSLWHLQGHÀQLGDmente. A pesar de que no existe una forma de determinar sen ξ(0.01), sabemos que todos los valores del seno se encuentran en el intervalo [−1, 1], por lo que el error que se presenta si utilizamos la aproximación 0.99995 para el valor de cos 0.01 está limitado por

| cos(0.01) − 0.99995| = 0.16 × 10−6 | sen ξ(0.01)| ≤ 0.16 × 10−6 . Por lo tanto, la aproximación 0.99995 corresponde por lo menos a los primeros cinco dígitos de cos 0.01 y

0.9999483 < 0.99995 − 1.6 × 10−6 ≤ cos 0.01 ≤ 0.99995 + 1.6 × 10−6 < 0.9999517. El límite del error es mucho más grande que el error real. Esto se debe, en parte, al escaso límite que usamos para | sen ξ(x)|. En el ejercicio 27 se muestra que para todos los valores de x, tenemos | sen x| ≤ |x|. Puesto que 0 ≤ ξ < 0.01, podríamos haber usado el hecho de que | sen ξ(x)| ≤ 0.01 en la fórmula de error, lo cual produce el límite 0.16 × 10−8 . b) Puesto que f (0) = 0, el tercer polinomio de Taylor con el término restante alrededor de x0 = 0 es

1 1 cos x = 1 − x 2 + x 4 cos ξ˜ (x), 2 24 donde 0 < ξ˜ (x) < 0.01. El polinomio de aproximación sigue siendo el mismo y la aproximación sigue siendo 0.99995, pero ahora tenemos mayor precisión. Puesto que | cos ξ˜ (x)| ≤ 1 para todas las x, obtenemos

1 4 1 x cos ξ˜ (x) ≤ (0.01)4 (1) ≈ 4.2 × 10−10 . 24 24 por lo tanto

| cos 0.01 − 0.99995| ≤ 4.2 × 10−10 , y

0.99994999958 = 0.99995 − 4.2 × 10−10 ≤ cos 0.01 ≤ 0.99995 + 4.2 × 10−10 = 0.99995000042. (OHMHPSORLOXVWUDORVGRVREMHWLYRVGHODQiOLVLVQXPpULFR i)

Encuentre una aproximación a la solución de un problema determinado.

ii)

Determine un límite o cota para la precisión de la aproximación.

/RVSROLQRPLRVGH7D\ORUHQDPEDVSDUWHVSURSRUFLRQDQODPLVPDUHVSXHVWDSDUDL SHURHO WHUFHURSURYHHXQDUHVSXHVWDPXFKRPHMRUSDUDLL TXHHOVHJXQGR7DPELpQSRGHPRVXWLOL]DU estos polinomios para obtener aproximaciones de las integrales.

1.2

Ilustración

Errores de redondeo y aritmética computacional

11

Podemos utilizar el tercer polinomio de Taylor y su término restante encontrado en el ejem0.1 SORSDUDDSUR[LPDU 0 cos x d x. Tenemos 0.1

0.1

cos x d x =

0

0

1 1 − x2 2 0.1

1 = x − x3 6

+ 0

dx + 1 24

1 1 = 0.1 − (0.1)3 + 6 24

1 24

0.1

0.1

x 4 cos ξ˜ (x) d x

0

x 4 cos ξ˜ (x) d x

0 0.1

x 4 cos ξ˜ (x) d x.

0

Por lo tanto, 0.1 0

1 cos x d x ≈ 0.1 − (0.1)3 = 0.09983. 6

Un límite o cota para el error en esta aproximación se determina a partir de la integral del término restante de Taylor y el hecho de que | cos ξ˜ (x)| ≤ 1 para todas las x:

1 24

0.1 0

1 x 4 cos ξ˜ (x) d x ≤ 24 ≤

1 24

0.1

x 4 | cos ξ˜ (x)| d x

0 0.1 0

x4 dx =

(0.1)5 = 8.3 × 10−8 . 120

El valor verdadero de esta integral es 0.1 0

0.1

cos x d x = sen x

= sen 0.1 ≈ 0.099833416647, 0

por lo que el error real para esta aproximación es 8.× 10−8, que se encuentra dentro del límite del error. La sección Conjunto de ejercicios 1.1 está disponible en línea. Encuentre la ruta de en las páginas preliminares.

1.2 Errores de redondeo y aritmética computacional La aritmética realizada con una calculadora o computadora es diferente a la aritmética que se imparte en los cursos √ de álgebra y cálculo. Podría esperarse que declaraciones como 2 +2 = 4, 4·8 = 32, y ( 3)2 = 3 siempre sean verdaderas; sin embargo con la aritmética computacional, esperamos resultados exactos para 2 + 2 = 4 y 4 · 8 = 32, pero no obtendre√ mos exactamente ( 3)2 = 3. Para comprender por qué esto es verdadero, debemos explorar HOPXQGRGHODDULWPpWLFDGHGtJLWRVÀQLWRV (QQXHVWURPXQGRPDWHPiWLFRWUDGLFLRQDOSHUPLWLPRVQ~PHURVFRQXQDFDQWLGDGLQÀQL√ ta de dígitos. La aritmética que usamos en este mundo GHÀQH 3 como el único número poVLWLYRTXHFXDQGRVHPXOWLSOLFDSRUVtPLVPRSURGXFHHOHQWHUR1RREVWDQWHHQHOPXQGR computacional, cada número representable sóORWLHQHXQQ~PHURÀMR\ÀQLWRGHGtJLWRV(VWR VLJQLÀFDTXHSRUHMHPSORVólo los números racionales, e incluso no todos ellos, se pueden √ representar de forma exacta. Ya que 3 no es racional, se proporciona una representación DSUR[LPDGDFX\RFXDGUDGRQRVHUiH[DFWDPHQWHDSHVDUGHTXHHVSUREDEOHTXHHVWpVXÀFLHQWHPHQWHFHUFDGHSDUDVHUDFHSWDEOHHQODPD\RUtDGHODVVLWXDFLRQHV(QWRQFHVHQPXchos casos, esta aritmética mecánica es satisfactoria y pasa sin importancia o preocupación, pero algunas veces surgen problemas debido a su discrepancia.

12

CAPÍTULO 1

Preliminares matemáticos y análisis de error

El error debido al redondeo debería esperarse siempre que se realizan cálculos con números que no son potencias de 2. Mantener este error bajo control es en extremo importante cuando el número de cálculos es grande.

El error que se produce cuando se utiliza una calculadora o computadora para realizar cálculos con números reales recibe el nombre de error de redondeo. Se presenta porque la DULWPpWLFDUHDOL]DGDHQXQDPiTXLQDLQFOX\HQ~PHURVFRQXQVRORQ~PHURÀQLWRGHGtJLWRV y esto da como resultado cálculos realizados únicamente con representaciones aproximadas de los números reales. En una computadora, sólo un subconjunto relativamente pequeño del sistema de números reales se usa para la representación de todos los números reales. Este subconjunto sólo contiene números racionales, tanto positivos como negativos, y almacena la parte fraccionaria, junto con una parte exponencial.

Números de máquina binarios En 1985, el Instituto para Ingenieros Eléctricos y Electrónicos (IEEE; Institute for ElectriFDODQG(OHFWURQLF(QJLQHHUV SXEOLFyXQUHSRUWHOODPDGRBinary Floating Point Arithmetic 6WDQGDUG ² (VWiQGDU SDUD OD DULWPpWLFD ELQDULD GH SXQWR ÁRWDQWH . En 2008 se publicó una versión actualizada con el nombre de IEEE 754-2008; la cual proporciona esWiQGDUHV SDUD Q~PHURV GH SXQWR ÁRWDQWH GHFLPDOHV \ ELQDULRV IRUPDWRV SDUD LQWHUFDPELR de datos, algoritmos para redondear operaciones aritméticas y manejo de excepciones. Se HVSHFLÀFDFXiOHVVRQORVIRUPDWRVSDUDODVSUHFLVLRQHVLQGLYLGXDOHVGREOHV\DPSOLDGDV\HQ JHQHUDOWRGRVORVIDEULFDQWHVGHPLFURFRPSXWDGRUDVTXHXWLOL]DQKDUGZDUHGHSXQWRÁRWDQWH siguen estos estándares. 8QDUHSUHVHQWDFLyQGHELWVGtJLWRELQDULR VHXVDSDUDXQQ~PHURUHDO(OSULPHUELW es un indicador de signo, denominado s. A éste le sigue un exponente de 11 bits, c, llamado característica, y una fracción binaria de 52 bits, f, llamada mantisa. La base para el exponente es 2. Puesto que los 52 dígitos binarios corresponden con dígitos decimales entre 16 y 17, podemos asumir que un número representado en este sistema tiene por lo menos 16 dígitos decimales de precisión. El exponente de 11 dígitos binarios provee un rango de 0 a 211 − 1 =6LQHPEDUJRXVDUVyORHQWHURVSRVLWLYRVSDUDHOH[SRQHQWHQRSHUPLWLUtDXQDUHpresentación adecuada de los números con magnitud pequeña. Para garantizar que estos Q~PHURVVRQLJXDOPHQWHUHSUHVHQWDEOHVVHUHVWDDODFDUDFWHUtVWLFDSRUORTXHHOUDQJR del exponente en realidad se encuentra entre −\ Para ahorrar almacenamiento y proporcionar una representación única para cada número GHSXQWRÁRWDQWHVHLPSRQHXQDQRUPDOL]DFLyQ3RUPHGLRGHHVWHVLVWHPDREWHQHPRVXQ Q~PHURGHSXQWRÁRWDQWHGHODIRUPD

(−1)s 2c−1023 (1 + f ). Ilustración

Considere el número de máquina 0 10000000011 1011100100010000000000000000000000000000000000000000. El bit más a la izquierda es s = 0, lo cual indica que es un número positivo. Los siguientes 11 bits, 10000000011, proveen la característica y son equivalentes al número decimal

c = 1 · 210 + 0 · 29 + · · · + 0 · 22 + 1 · 21 + 1 · 20 = 1024 + 2 + 1 = 1027. La parte exponencial del número es, por lo tanto, 21027−1023 = 24/RVELWVÀQDOHVHVSHFLÀFDQTXHODPDQWLVDHV

f =1·

1 2

1

+1·

1 2

3

+1·

1 2

4

+1·

1 2

5

+1·

1 2

8

+1·

1 2

12

.

1.2

Errores de redondeo y aritmética computacional

13

Como secuencia, este número de máquina representa precisamente el número decimal

(−1)s 2c−1023 (1 + f ) = (−1)0 · 21027−1023 1 +

1 1 1 1 1 1 + + + + + 2 8 16 32 256 4096

= 27.56640625. Sin embargo, el siguiente número de máquina más pequeño es 0 10000000011 1011100100001111111111111111111111111111111111111111, el siguiente número de máquina más grande es 0 10000000011 1011100100010000000000000000000000000000000000000001. (VWRVLJQLÀFDTXHQXHVWURQ~PHURGHPiTXLQDRULJLQDOQRVóORUHSUHVHQWDVLQR WDPELpQODPLWDGGHORVQ~PHURVUHDOHVTXHVHHQFXHQWUDQHQWUH\HOVLJXLHQWH Q~PHURGHPiTXLQDPiVSHTXHxRDVtFRPRODPLWDGGHORVQ~PHURVHQWUH\ el siguiente número de máquina más grande. Para ser preciso, representa cualquier número > El número positivo normalizado más pequeño que se puede representar tiene s = 0, c = 1, y f = 0 y es equivalente

2−1022 · (1 + 0) ≈ 0.22251 × 10−307 , y el más grande tiene s = 0, c =\f = 1 − 2−52 y es equivalente a

21023 · (2 − 2−52 ) ≈ 0.17977 × 10309 . Los números que se presentan en los cálculos que tienen una magnitud menor que

2−1022 · (1 + 0) resultan en un subdesbordamiento\HQJHQHUDOVHFRQÀJXUDQHQFHUR/RVQ~PHURVVXSHriores a

21023 · (2 − 2−52 ) resultan en desbordamiento y, comúnmente, causan que los cálculos se detengan (a menos TXH HO SURJUDPD KD\D VLGR GLVHxDGR SDUD GHWHFWDU HVWDV SUHVHQFLDV 2EVHUYH TXH H[LVWHQ dos representaciones para el número cero: un 0 positivo cuando s = 0, c = 0 y f = 0, y un 0 negativo cuando s = 1, c = 0 y f = 0.

Números de máquina decimales (OXVRGHGtJLWRVELQDULRVWLHQGHDRFXOWDUODVGLÀFXOWDGHVFRPSXWDFLRQDOHVTXHVHSUHVHQWDQ FXDQGRVHXVDXQDFROHFFLyQÀQLWDGHQ~PHURVGHPiTXLQDSDUDUHSUHVHQWDUWRGRVORVQ~PHros reales. Para examinar estos problemas, utilizaremos números decimales más familiares HQOXJDUGHXQDUHSUHVHQWDFLyQELQDULD(QHVSHFtÀFRVXSRQHPRVTXHORVQ~PHURVPiTXLQD VHUHSUHVHQWDQHQIRUPDWRQRUPDOL]DGRGHSXQWRÁRWDQWHdecimal

±0.d1 d2 . . . dk × 10n ,

1 ≤ d1 ≤ 9,

y

0 ≤ di ≤ 9,

14

CAPÍTULO 1

Preliminares matemáticos y análisis de error

para cada i = 2, . . . , k. Los números de esta forma reciben el nombre de números de máquina decimales de dígito k. Cualquier número real positivo dentro del rango numérico de la máquina puede ser normalizado a la forma

y = 0.d1 d2 . . . dk dk+1 dk+2 . . . × 10n . El error que resulta de reemplazar un número con HVWDIRUPDGHSXQWRÁRWDQWH se llama error de redondeo, independientemente de si se usa el método de redondeo o de corte.

/DIRUPDGHSXQWRÁRWDQWHGH\TXHVHGHQRWDfl(y , se obtiene al terminar la mantisa de y en los dígitos decimales de k. Existen dos maneras comunes para realizar esta terminación. Un método, llamado de corte, es simplemente cortar los dígitos dk+1 dk+2 . . . Esto produce ODIRUPDGHSXQWRÁRWDQWH

f l(y) = 0.d1 d2 . . . dk × 10n . El otro método, llamado redondeo, suma 5 × 10n−(k+1) a y y entonces corta el resultado para obtener un número con la forma

f l(y) = 0.δ1 δ2 . . . δk × 10n . Para redondear, cuando dk+1 ≥ 5, sumamos 1 a dk para obtener fl(y HVGHFLUredondeamos hacia arriba. Cuando dk+1 < 5, simplemente cortamos todo, excepto los primeros dígitos k; es decir, redondeamos hacia abajo. Si redondeamos hacia abajo, entonces δi = di, para cada i = 1, 2, . . . , k. Sin embargo, si redondeamos hacia arriba, los dígitos (e incluso el H[SRQHQWH SXHGHQFDPELDU Ejemplo 1

'HWHUPLQHORVYDORUHVD GHFRUWH\E GHUHGRQGHRGHFLQFRGtJLWRVGHOQ~PHURLUUDFLRQDOπ. Solución El número πWLHQHXQDH[SDQVLyQGHFLPDOLQÀQLWDGHODIRUPDπ =.... Escrito en una forma decimal normalizada, tenemos

π = 0.314159265 . . . × 101 . En general, el error relativo es una mejor medición de precisión que el error absoluto porque considera el tamaño del número que se va a aproximar.

a) (OIRUPDWRGHSXQWRÁRWDQWHGHπ usando el recorte de cinco dígitos es

f l(π ) = 0.31415 × 101 = 3.1415. b) El sexto dígito de la expansión decimal de π es un 9, por lo que el formato de punto ÁRWDQWHGHπ con redondeo de cinco dígitos es

f l(π ) = (0.31415 + 0.00001) × 101 = 3.1416. /DVLJXLHQWHGHÀQLFLyQGHVFULEHWUHVPpWRGRVSDUDPHGLUHUURUHVGHDSUR[LPDFLyQ Definición 1.15

Suponga que p ∗ es una aproximación a p. El error real es p − p ∗, el error absoluto es | p − p∗ | | p − p ∗ |, y el error relativo es , siempre y cuando p = 0. | p| Considere los errores real, absoluto y relativo al representar p con p ∗ en el siguiente ejemplo.

Ejemplo 2

Determine los errores real, absoluto y relativo al aproximar p con p ∗ cuando

a)

p = 0.3000 × 101 y p ∗ = 0.3100 × 101 ;

b)

p = 0.3000 × 10−3 y p ∗ = 0.3100 × 10−3 ;

c)

p = 0.3000 × 104 y p ∗ = 0.3100 × 104 .

1.2 Solución

a) b) A menudo no podemos encontrar un valor preciso para el error verdadero en una aproximación. Por el contrario, encontramos una cota para el error, lo cual nos proporciona un error del “peor caso”.

c)

15

Errores de redondeo y aritmética computacional

Para p =× 101 y p ∗ =× 101, el error real es <0.1, el error absoluto es 0.1 y el error relativo es 0.3333 × 10−1. Para p =× 10ï y p ∗ =× 10ï, el error real es <0.1 × 10ï, el error absoluto es 0.1 × 10ï y el error relativo es 0.3333 × 10−1. Para p =× 10 y p ∗ =×HOHUURUUHDOHVï× 10, el error absoluto es 0.1 × 10 y, de nuevo, el error relativo es 0.3333 × 10−1 .

Este ejemplo muestra que el mismo error relativo, 0.3333 × 10−1, se presenta para errores absolutos ampliamente variables. Como una medida de precisión, el error absoluto puede ser HQJDxRVR\HOHUURUUHODWLYRPiVVLJQLÀFDWLYRGHELGRDTXHHVWHHUURUFRQVLGHUDHOWDPDxR del valor. Un límite de error es un número no negativo mayor que el error absoluto. Algunas veces se obtiene con los métodos de cálculo para encontrar el valor absoluto máximo de una IXQFLyQ(VSHUDPRVHQFRQWUDUHOOtPLWHVXSHULRUPiVSHTXHxRSRVLEOHSDUDHOHUURUDÀQGH obtener un estimado del error real que es lo más preciso posible. /D VLJXLHQWH GHÀQLFLyQ XVD HO HUURU UHODWLYR SDUD SURSRUFLRQDU XQD PHGLGD GH GtJLWRV VLJQLÀFDWLYRVGHSUHFLVLyQSDUDXQDDSUR[LPDFLyQ

Definición 1.16 A menudo, el término dígitos VLJQLÀFDWLYRV se usa para describir vagamente el número de dígitos decimales que parecen VHUH[DFWRV/DGHÀQLFLyQHVPiV precisa y provee un concepto continuo.

Tabla 1.1

Se dice que el número p ∗ se aproxima a p para t GtJLWRVVLJQLÀFDWLYRVRFLIUDV VLt es el entero no negativo más grande para el que

| p − p∗ | ≤ 5 × 10−t . | p| /DWDEODLOXVWUDODQDWXUDOH]DFRQWLQXDGHORVGtJLWRVVLJQLÀFDWLYRVDOHQXPHUDUSDUD los diferentes valores de p, el límite superior mínimo de | p − p ∗ |, denominado máx. | p − p ∗ |, cuando p ∗ concuerda con pHQFXDWURGtJLWRVVLJQLÀFDWLYRV

p

0.1

0.5

100

1000

5000

9990

10000

máx | p − p ∗ |

0.00005

0.00025

0.05

0.5

2.5

4.995

5.

Al regresar a la representación de los números de máquina, observamos que la represenWDFLyQGHSXQWRÁRWDQWHf l(y SDUDHOQ~PHURy tiene el error relativo

y − f l(y) . y Si se usan k dígitos decimales y corte para la representación de máquina de

y = 0.d1 d2 . . . dk dk+1 . . . × 10n , entonces

y − f l(y) 0.d1 d2 . . . dk dk+1 . . . × 10n − 0.d1 d2 . . . dk × 10n = y 0.d1 d2 . . . × 10n =

0.dk+1 dk+2 . . . 0.dk+1 dk+2 . . . × 10n−k = × 10−k . n 0.d1 d2 . . . × 10 0.d1 d2 . . .

CAPÍTULO

Interpolación y aproximación polinomial

Introducción Se realiza un censo de la población de Estados Unidos cada 10 años. La siguiente tabla muestra la población, en miles de personas, desde 1960 hasta 2010, y los datos también se UHSUHVHQWDQHQODÀJXUD Año Población (en miles)

1960

1970

1980

1990

2000

2010

179 323

203 302

226542

249 633

281 422

308 746

P(t) 3 3 10 8

2 3 10 8 Población

3

1 3 10 8

1960 1970 1980 1990 2000 2010 Año

t

Al revisar estos datos, podríamos preguntar si se podrían usar para efectuar un cálculo razonable de la población, digamos, en 1975 o incluso en el año 2020. Las predicciones de este tipo pueden obtenerse por medio de una función que se ajuste a los datos proporcionados. Este proceso recibe el nombre de interpolación y es el tema de este capítulo. Este problema de población se considera a lo largo del capítulo y en los ejercicios 19 de la sección 3.1, 17 de la sección 3.3 y 24 de la sección 3.5.

77

78

CAPÍTULO 3

Interpolación y aproximación polinomial

3.1 Interpolación y el polinomio de Lagrange Una de las clases más útiles y conocidas de funciones que mapean el conjunto de números reales en sí mismo son los polinomios algebraicos, el conjunto de funciones de la forma

Pn (x) = an x n + an−1 x n−1 + · · · + a1 x + a0 , donde n es un entero positivo y a0, 7, an son constantes reales. Una razón de su importancia es que se aproximan de manera uniforme a las funciones continuas. Con esto queremos decir TXHGDGDXQDIXQFLyQGHÀQLGD\FRQWLQXDVREUHXQLQWHUYDORFHUUDGR\DFRWDGRH[LVWHXQ polinomio que está tan “cerca” de la función dada como se desee. Este resultado se expresa FRQSUHFLVLyQHQHOWHRUHPDGHDSUR[LPDFLyQGH:HLHUVWUDVVFRQVXOWHODÀJXUD

Figura 3.1 y

y 5 f (x) 1 y 5 P (x) y 5 f (x) y 5 f (x) 2

a

Teorema 3.1

b

x

(Teorema de aproximación de Weierstrass) Suponga que fHVWiGHÀQLGD\HVFRQWLQXDHQ>a, b]. Para cada . 0, existe un polinomio P (x), con la propiedad de que

A menudo, se hace referencia a Karl Weierstrass (1815–1897) como el padre del análisis moderno debido a su insistencia sobre el rigor en la demostración de resultados matemáticos. Fue fundamental para el desarrollo de pruebas de convergencia de series y para determinar formas GHGHÀQLUULJXURVDPHQWHORV números irracionales. Fue el primero en demostrar que una función podría ser continua en todas partes, pero diferenciable en ninguna parte, un resultado que escandalizó a algunos de sus contemporáneos.

| f (x) − P(x)|

para todas las x en [a, b].

La prueba de este teorema se puede encontrar en la mayoría de los textos básicos sobre DQiOLVLVUHDOFRQVXOWHSRUHMHPSOR>%DUW@SS² Otra razón importante para considerar la clase de polinomios en la aproximación de IXQFLRQHVHVTXHODGHULYDGD\ODLQWHJUDOLQGHÀQLGDGHXQSROLQRPLRVRQIiFLOHVGHGHWHUPLnar y también son polinomios. Por esta razón, a menudo se usan polinomios para aproximar funciones continuas. Los polinomios de Taylor se presentaron en la sección 1.1, donde se describieron como uno de los componentes básicos del análisis numérico. Debido a su importancia, se podría esperar que la aproximación polinomial usará estas funciones en gran medida; sin embargo, éste no es el caso. Los polinomios de Taylor concuerdan tanto como es posible con una IXQFLyQGDGDHQXQSXQWRHVSHFtÀFRSHURFRQFHQWUDQVXSUHFLVLyQFHUFDGHHVHSXQWR8Q buen polinomio de aproximación debe dar precisión relativa sobre un intervalo completo y, en general, los polinomios de Taylor no lo hacen. Por ejemplo, suponga que calculamos los

3.1

Interpolación y el polinomio de Lagrange

79

primeros seis polinomios de Taylor alrededor de x0 5 0 para f (x) 5 e x. Ya que las derivadas de f (x) son todas ex, que evaluadas en x0 5 0 dan 1, los polinomios de Taylor son Se publicó muy poco del trabajo de Weierstrass durante su vida; no obstante, sus conferencias, en especial sobre la teoría de las IXQFLRQHVLQÁX\HURQGHPDQHUD VLJQLÀFDWLYDHQXQDJHQHUDFLyQ completa de estudiantes.

P0 (x) = 1,

P1 (x) = 1 + x,

P4 (x) = 1 + x +

x2 , 2

P2 (x) = 1 + x +

x3 x4 x2 + + , 2 6 24

P3 (x) = 1 + x +

P5 (x) = 1 + x +

y

x3 x2 + , 2 6

x3 x4 x5 x2 + + + . 2 6 24 120

/DVJUiÀFDVGHORVSROLQRPLRVVHPXHVWUDQHQODÀJXUDREVHUYHTXHLQFOXVRSDUD los polinomios de grado más alto, el error empeora progresivamente conforme nos alejamos de cero).

Figura 3.2 y 20

y 5 P5(x)

y 5 ex

y 5 P4(x)

15

y 5 P3(x) 10

y 5 P2(x)

5

y 5 P1(x) y 5 P0(x)

21

2

1

x

3

Aunque se obtienen mejores aproximaciones para f (x) 5 e x si se usan polinomios de Taylor, esto no es verdad para todas las funciones. Considere, como un ejemplo extremo, usar la expansión en polinomios de Taylor de diferentes grados para f (x) 5 1/ x alrededor de x0 5 1 para aproximar f (3) 5 1/3. Puesto que

f (x) = x −1 , f (x) = −x −2 , f (x) = (−1)2 2 · x −3 , y, en general,

f (k) (x) = (−1)k k!x −k−1 , los polinomios de Taylor son n

Pn (x) = k=0

f (k) (1) (x − 1)k = k!

n

(−1)k (x − 1)k .

k=0

Para aproximar f (3) 5 1/3 mediante P n (3) para valores cada vez mayores de n, obtenemos los valores en la tabla 3.1 (¡un terrible fracaso!). Cuando aproximamos f (3) 5 1/3 mediante P n (3) y para valores más grandes de n, la aproximación se vuelve cada vez más imprecisa.

Tabla 3.1

n

0

1

2

3

4

5

6

7

Pn (3)

1

−1

3

−5

11

−21

43

−85

80

CAPÍTULO 3

Interpolación y aproximación polinomial

Para los polinomios de Taylor, toda la información que se usa en la aproximación se concentra en el único número x0, por lo que, en general, éstos darán aproximaciones imprecisas conforme nos alejamos de x0. Esto limita la aproximación de polinomios de Taylor a situaciones en las que las aproximaciones sólo se necesitan en números cercanos a x0. Para SURSyVLWRVFRPSXWDFLRQDOHVRUGLQDULRVHVPiVHÀFLHQWHXVDUPpWRGRVTXHLQFOX\DQLQIRUmación en varios puntos. Consideramos esto en el resto del capítulo. El uso principal de los polinomios de Taylor en el análisis numérico no tiene propósitos de aproximación, sino la derivación de técnicas numéricas y el cálculo de errores.

Polinomios de interpolación de Lagrange El problema de determinar un polinomio de grado uno que pasa por diferentes puntos (x0, y0) y (x1, y1) es igual al de aproximar una función f para la que f(x0 ) 5 y0 y f (x1 ) 5 y1 por medio de un polinomio de primer grado que se interpola, o que coincida con los valores de f en los puntos determinados. El uso de estos polinomios para aproximación dentro del intervalo GHWHUPLQDGRPHGLDQWHSXQWRVÀQDOHVUHFLEHHOQRPEUHGHinterpolación. 'HÀQDODVIXQFLRQHV

L 0 (x) =

x − x1 x0 − x1

y L 1 (x) =

x − x0 . x1 − x0

El polinomio de interpolación de Lagrange lineal a través de (x0, y0) y (x1, y1) es

P(x) = L 0 (x) f (x0 ) + L 1 (x) f (x1 ) =

x − x1 x − x0 f (x0 ) + f (x1 ). x0 − x1 x1 − x0

Observe que

L 0 (x0 ) = 1,

L 0 (x1 ) = 0,

L 1 (x0 ) = 0,

y L 1 (x1 ) = 1,

lo cual implica que

P(x0 ) = 1 · f (x0 ) + 0 · f (x1 ) = f (x0 ) = y0 y

P(x1 ) = 0 · f (x0 ) + 1 · f (x1 ) = f (x1 ) = y1 . Por lo que P es el único polinomio de grado a lo más 1 que pasa por (x0, y0) y (x1, y1). Ejemplo 1

Determine el polinomio de interpolación de Lagrange que pasa por los puntos (2, 4) y (5, 1). Solución en este caso, tenemos

L 0 (x) =

1 x −5 = − (x − 5) y 2−5 3

L 1 (x) =

1 x −2 = (x − 2), 5−2 3

por lo que 1 4 20 1 2 1 + x − = −x + 6. P(x) = − (x − 5) · 4 + (x − 2) · 1 = − x + 3 3 3 3 3 3 /DJUiÀFDGHy 5 P(x VHPXHVWUDHQODÀJXUD

3.1

Interpolación y el polinomio de Lagrange

81

Figura 3.3 y (2,4)

4 3 2

y 5 P(x) = 2x 1 6

1

1

2

3

4

(5,1)

5

x

Para generalizar el concepto de interpolación lineal, considere la construcción de un polinomio de grado n que pasa a través de n 1 1 puntos

(x0 , f (x0 )), (x1 , f (x1 )), . . . , (xn , f (xn )). 9pDVHODÀJXUD

Figura 3.4 y

y 5 f (x) y 5 P(x)

x0

x1

x2

xn

x

En este caso, primero construimos, para cada k 5 0, 1, 7, n, una función L n,k (x) con la propiedad de que Ln,k (xi ) = 0 cuando i = k y L n,k (x k) 5 1. Para satisfacer Ln,k (xi ) = 0 para cada i = k se requiere que el numerador de L n,k (x) contenga el término

(x − x0 )(x − x1 ) · · · (x − xk−1 )(x − xk+1 ) · · · (x − xn ). Para satisfacer L n,k (x k) 5 1, el denominador de L n,k (x) debe ser el mismo término, pero evaluado en x 5 xk. Por lo tanto,

L n,k (x) =

(x − x0 ) · · · (x − xk−1 )(x − xk+1 ) · · · (x − xn ) . (xk − x0 ) · · · (xk − xk−1 )(xk − xk+1 ) · · · (xk − xn )

8QERVTXHMRGHODJUiÀFDGHXQDL n,k (cuando n HVSDU VHPXHVWUDHQODÀJXUD

82

CAPÍTULO 3

Interpolación y aproximación polinomial

Figura 3.5 L n,k(x)

1

x0

x1

...

x k21

xk

x k11

...

x n21

xn

x

El polinomio de interpolación se describe fácilmente una vez que se conoce la forma L n,k . Este polinomio, llamado enésimo polinomio de interpolación de Lagrange,VHGHÀQH en el siguiente teorema. Teorema 3.2 La fórmula de interpolación nombrada por Joseph Louis Lagrange (1736–1813) probablemente era conocida por Newton alrededor de 1675, pero al parecer fue publicada por primera vez en 1779 por Edward Waring (1736–1798). Lagrange escribió mucho sobre el tema de interpolación y su trabajo tuvo XQDLQÁXHQFLDVLJQLÀFDWLYDVREUH los matemáticos posteriores. Él publicó este resultado en 1795.

Si x0, x1, 7, x n son n 1 1 números distintos y f es una función cuyos valores están determinados en estos números, entonces existe un único polinomio P(x) de grado a lo sumo n con

f (xk ) = P(xk ), Este polinomio está determinado por

n

P(x) = f (x0 )L n,0 (x) + · · · + f (xn )L n,n (x) =

= a1 ∗ a2 ∗ a3 .

f (xk )L n,k (x),

(3.1)

k=0

donde, para cada k = 0, 1, . . . , n, L n,k (x) =

El símbolo se usa para escribir productos de manera compacta y es similar al símbolo , que se utiliza para escribir sumas. Por ejemplo 3 i=0 ai

para cada k = 0, 1, . . . , n.

(x − x0 )(x − x1 ) · · · (x − xk−1 )(x − xk+1 ) · · · (x − xn ) (xk − x0 )(xk − x1 ) · · · (xk − xk−1 )(xk − xk+1 ) · · · (xk − xn ) n

= i=0 i =k

(3.2)

(x − xi ) . (xk − xi )

Escribiremos Ln,k (x) simplemente como L k (x) cuando no haya confusión en cuanto a su grado.

Ejemplo 2

a)

Use los números (llamados nodos) x0 5 2, x1 5 2.75 y x2 5 4 para encontrar el polinomio de interpolación de Lagrange de segundo grado para f (x) 5 1/x.

b)

Use este polinomio para aproximar f(3) 5 1/3.

Solución a)3ULPHURGHWHUPLQDPRVORVFRHÀFLHQWHVSROLQyPLFRVL 0 (x), L 1 (x) y L 2 (x). En forma anidada, estos son

L 0 (x) =

(x − 2.75)(x − 4) 2 = (x − 2.75)(x − 4), (2 − 2.75)(2 − 4) 3

L 1 (x) =

16 (x − 2)(x − 4) = − (x − 2)(x − 4), (2.75 − 2)(2.75 − 4) 15

L 2 (x) =

2 (x − 2)(x − 2.75) = (x − 2)(x − 2.75). (4 − 2)(4 − 2.75) 5

y

3.1

Interpolación y el polinomio de Lagrange

83

Además, f (x0 ) = f (2) = 1/2, f (x1 ) = f (2.75) = 4/11, y f (x2 ) = f (4) = 1/4, por lo que 2

P(x) =

f (xk )L k (x) k=0

1 64 1 (x − 2.75)(x − 4) − (x − 2)(x − 4) + (x − 2)(x − 2.75) 3 165 10 49 1 2 35 x − x+ . = 22 88 44 =

b) Una aproximación para f (3) = 1/3 (véase la figura 3.6) es f (3) ≈ P(3) =

105 49 29 9 − + = ≈ 0.32955. 22 88 44 88

Recuerde que en la sección de apertura de este capítulo (consulte la tabla 3.1), encontramos que ninguna expansión en polinomios de Taylor alrededor de x0 5 1 se puede usar para aproximar razonablemente f(x) 5 1/x en x 5 3. Figura 3.6 y 4 3 2

y 5 f (x)

1 y 5 P(x) 1

2

3

4

5

x

El siguiente paso es calcular un residuo o cota para el error involucrado en la aproximación de una función mediante un polinomio de interpolación. Teorema 3.3

Existen otras formas de expresar el término de error para el polinomio de Lagrange, pero ésta puede ser la forma más útil y la que concuerda más estrechamente con la forma de error del polinomio estándar de Taylor.

Suponga x0 , x1 , . . . , xnVRQQ~PHURVGLVWLQWRVHQHOLQWHUYDOR>a, b] y f ∈ C n+1 [a, b]. Entonces, para cada x en >a, b], existe un número ξ(x) (generalmente no conocido) entre mín {x0, x1, 7, xn} y máx{x0, x1, 7, xn} y, por lo tanto, en (a, b), con

f (x) = P(x) +

f (n+1) (ξ(x)) (x − x0 )(x − x1 ) · · · (x − xn ), (n + 1)!

(3.3)

donde P(x) es el polinomio de interpolación determinado en la ecuación (3.1). Demostración Primero observe que si x 5 x k para cualquier k 5 0, 1, 7, n, entonces f(xk ) 5 P(x k ) y al elegir ξ (x k )de manera arbitraria en (a, b) se obtiene la ecuación (3.3).

84

CAPÍTULO 3

Interpolación y aproximación polinomial

Si x = xk , para todas las k = 0, 1, . . . , nGHÀQDODIXQFLyQg para tHQ>a, b] mediante

g(t) = f (t) − P(t) − [ f (x) − P(x)]

(t − x0 )(t − x1 ) · · · (t − xn ) (x − x0 )(x − x1 ) · · · (x − xn ) n

= f (t) − P(t) − [ f (x) − P(x)] i=0

(t − xi ) . (x − xi )

Puesto que f ∈ C n+1 [a, b], y P ∈ C ∞ [a, b], se sigue que g ∈ C n+1 [a, b]. Para t = xk , tenemos n

g(xk ) = f (xk ) − P(xk ) − [ f (x) − P(x)] i=0

(xk − xi ) = 0 − [ f (x) − P(x)] · 0 = 0. (x − xi )

Además, n

g(x) = f (x) − P(x) − [ f (x) − P(x)] i=0

(x − xi ) = f (x) − P(x) − [ f (x) − P(x)] = 0. (x − xi )

Por lo tanto, g ∈ C n+1 [a, b], y g se anula en los n 1 2 números distintos x, x0 , x1 , . . . , xn. Por el teorema generalizado de Rolle 1.10, existe un número ξ en (a, b) para el que g (n+1) (ξ ) = 0.Por lo que,

0 = g (n + 1) (ξ ) = f (n+1) (ξ ) − P (n+1) (ξ ) − [ f (x) − P(x)]

d n+1 dt n+1

n i=0

(t − xi ) (x − xi )

. (3.4) t=ξ

Sin embargo, P(x) es un polinomio de grado a lo sumo n, por lo que la derivada (n 1 1), n P (n+1) (x), es cero. Además i=0 [(t − xi )/(x − xi )] es un polinomio de grado (n 1 1), por lo que n i=0

(t − xi ) = (x − xi )

1 n i=0 (x

− xi )

t n+1 + (términos de menor grado en t),

y d n+1 dt n+1

n i=0

(t − xi ) = (x − xi )

(n + 1)! . − xi )

n i=0 (x

Ahora, la ecuación (3.4) se convierte en 0 = f (n+1) (ξ ) − 0 − [ f (x) − P(x)]

(n + 1)! , − xi )

n i=0 (x

y, después de resolver f (x), tenemos f (x) = P(x) +

f (n+1) (ξ ) (n + 1)!

n

(x − xi ). i=0

La fórmula de error en el teorema 3.3 es un resultado teórico importante porque los polinomios de Lagrange se usan ampliamente para deducir la diferenciación numérica y los métodos de integración. Las cotas de error para estas técnicas se obtienen a partir de la fórmula del error de Lagrange. Observe que la forma del error para el polinomio de Lagrange es bastante similar a la del polinomio de Taylor. El enésimo polinomio de Taylor alrededor de x0 concentra toda la información conocida en x0 y tiene un término de error de la forma

f (n+1) (ξ(x)) (x − x0 )n+1 . (n + 1)!

3.1

Interpolación y el polinomio de Lagrange

85

El polinomio de Lagrange de grado n utiliza información en los distintos números x0, x1, , xn y, en lugar de (x 2 x0)n su fórmula de error utiliza el producto de los n 1 1 términos (x − x0 ), (x − x1 ), . . . , (x − xn ):

f (n+1) (ξ(x)) (x − x0 )(x − x1 ) · · · (x − xn ). (n + 1)! Ejemplo 3

En el ejemplo 2 encontramos el segundo polinomio de Lagrange para f(x) 5 1/x HQ>@ usando los nodos x0 5 2, x1 5 2.75 y x2 = 4. Determine la forma del error para este polinomio y el error máximo cuando el polinomio se usa para aproximar f (x) para x ∈ >@ Como f (x) = x −1, tenemos

Solución

f (x) = −x −2 ,

f (x) = 2x −3 ,

y

f (x) = −6x −4 .

En consecuencia, el segundo polinomio de Lagrange tiene el error de la forma

f (ξ(x)) (x− x0 )(x− x1 )(x− x2 ) = − (ξ(x))−4 (x− 2)(x− 2.75)(x− 4), para ξ(x)en(2, 4). 3! El valor máximo de (ξ(x))−4 en el intervalo es 2−4 = 1/16. Ahora necesitamos determinar el valor máximo en este intervalo del valor absoluto del polinomio

g(x) = (x − 2)(x − 2.75)(x − 4) = x 3 −

35 2 49 x + x − 22. 4 2

Como Dx

x3 −

35 2 49 x + x − 22 4 2

= 3x 2 −

49 1 35 x+ = (3x − 7)(2x − 7), 2 2 2

los puntos críticos se presentan en x=

7 , con g 3

7 3

=

25 , 108

y

x=

7 , con g 2

7 2

=−

9 . 16

Por lo tanto, el error máximo es 1 9 9 f (ξ(x)) |(x − x0 )(x − x1 )(x − x2 )| ≤ ≈ 0.03515625. − = 3! 16 16 256 El siguiente ejemplo ilustra cómo se puede usar la fórmula del error para preparar una tabla de datos que garantizará un error de interpolación dentro de una cota establecida. Ejemplo 4

Suponga que se va a preparar una tabla para la función f (x) = e x , para x en [0, 1]. Imagine que el número de lugares decimales proporcionado por entrada es d $ 8 y que h, el tamaño del paso es la diferencia entre valores adyacentes x. ¿Qué tamaño de paso h garantizará que la interpolación lineal proporcione un error absoluto a lo máximo de 10–6 para todas las x HQ>@" Solución

Sean x0, x1, 7 los números en los que se evalúa f y xHVWiHQ>@\VXSRQJDTXHj satisface xj # x # x j 11. La ecuación (3.3) implica que el error en la interpolación lineal es

| f (x) − P(x)| =

f (2) (ξ ) | f (2) (ξ )| (x − x j )(x − x j+1 ) = |(x − x j )||(x − x j+1 )|. 2! 2

Como el tamaño del paso es h, entonces x j = j h, x j+1 = ( j + 1)h, y | f (x) − P(x)| ≤

| f (2) (ξ )| |(x − j h)(x − ( j + 1)h)|. 2!

86

CAPÍTULO 3

Interpolación y aproximación polinomial

Por lo tanto,

máxξ ∈[0,1] eξ máx |(x − j h)(x − ( j + 1)h)| x j ≤x≤x j+1 2 e ≤ máx |(x − j h)(x − ( j + 1)h)|. 2 x j ≤x≤x j+1

| f (x) − P(x)| ≤

Considere la función g(x) = (x − j h)(x − ( j + 1)h), para j h ≤ x ≤ ( j + 1)h. Luego g (x) = (x − ( j + 1)h) + (x − j h) = 2 x − j h −

h 2

,

el único punto crítico para g se encuentra en x = j h + h/2, con g( j h + h/2) = (h/2)2 = h 2 /4. Puesto que g( j h) = 0 y g(( j + 1)h) = 0, el valor máximo de |g (x)| en [j h, ( j + 1)h] se debe presentar en el punto crítico, lo cual implica que (véase el ejercicio 21)

| f (x) − P(x)| ≤

e 2

máx

x j ≤x≤x j+1

|g(x)| ≤

e h2 eh 2 · = . 2 4 8

Por consiguiente, para garantizar que el error en la interpolación lineal está acotado por 1026HVVXÀFLHQWHHOHJLUh de tal forma que

eh 2 ≤ 10−6 . 8

Esto implica que

h < 1.72 × 10−3 .

Puesto que n 5 (1 2 0)/h debe ser un entero, una selección razonable para el tamaño del paso es h 5 0.001. La sección Conjunto de ejercicios 3.1 está disponible en línea. Encuentre la ruta de en las páginas preliminares.

3.2 Aproximación de datos y método de Neville En la sección anterior encontramos una representación explícita para los polinomios de Lagrange y su error cuando se aproxima una función sobre un intervalo. El uso frecuente de estos polinomios implica la interpolación de datos tabulados. En este caso, una representación explícita del polinomio podría no ser necesaria, sólo los valores del polinomio en SXQWRVHVSHFtÀFRV(QHVWDVLWXDFLyQVHUtDSRVLEOHTXHODIXQFLyQVXE\DFHQWHDORVGDWRVQR se conozca, por lo que la forma explícita del error no se puede usar. Ahora, ilustraremos una aplicación práctica de interpolación en dicha situación. Ilustración

Tabla 3.2 x

f (x)

1.0 1.3 1.6 1.9 2.2

0.7651977 0.6200860 0.4554022 0.2818186 0.1103623

La tabla 3.2 lista los valores de una función f en diferentes puntos. Las aproximaciones para f (1.5) obtenidas con distintos polinomios de Lagrange que usan estos datos se comparará para probar y determinar la precisión de la aproximación. El polinomio lineal más apropiado usa x0 5 1.3 y x1 5 1.6 porque 1.5 se encuentra entre 1.3 y 1.6. El valor del polinomio de interpolación en 1.5 es

P1 (1.5) = =

(1.5 − 1.3) (1.5 − 1.6) f (1.3) + f (1.6) (1.3 − 1.6) (1.6 − 1.3) (1.5 − 1.3) (1.5 − 1.6) (0.6200860) + (0.4554022) = 0.5102968. (1.3 − 1.6) (1.6 − 1.3)

3.2

Aproximación de datos y método de Neville

87

Es posible usar razonablemente dos polinomios de grado dos, uno con x0 5 1.3, x1 5 1.6 y x2 5 1.9, lo cual nos da

P2 (1.5) =

(1.5 − 1.3)(1.5 − 1.9) (1.5 − 1.6)(1.5 − 1.9) (0.6200860) + (0.4554022) (1.3 − 1.6)(1.3 − 1.9) (1.6 − 1.3)(1.6 − 1.9) +

(1.5 − 1.3)(1.5 − 1.6) (0.2818186) = 0.5112857, (1.9 − 1.3)(1.9 − 1.6)

y uno con x0 5 1.0, x1 5 1.3 y x2 5 1.6, lo cual nos da Pˆ 2 (1.5) = 0.5124715. En el caso de tercer grado, también hay dos opciones razonables para el polinomio, una con x0 5 1.3, x1 5 1.6, x2 5 1.9 y x3 5 2.2, lo cual nos da P3(1.5) 5 0.5118302. La segunda aproximación de tercer grado se obtiene con x0 51.0, x1 5 1.3, x2 5 1.6 y x3 5 1.9, lo cual nos da Pˆ 3 (1.5) = 0.5118127. El polinomio de Lagrange de cuarto grado usa todas las entradas en la tabla. Con x0 5 1.0, x1 5 1.3, x2 = 1.6, x3 5 1.9 y x4 5 2.2, la aproximación es P4(1.5) = 0.5118200. Puesto que P3(1.5), Pˆ 3 (1.5) y P4(1.5) concuerdan con una exactitud de 2 3 1025 unidades, esperamos este grado de precisión para estas aproximaciones. También esperamos que P4(1.5) sea la aproximación más precisa ya que usa la mayor parte de los datos proporcionados. /DIXQFLyQTXHHVWDPRVDSUR[LPDQGRHVHQUHDOLGDGODIXQFLyQGH%HVVHOGHSULPHUD clase de orden cero, cuyo valor en 1.5 se conoce como 0.5118277. Por lo tanto, las verdaderas precisiones de las aproximaciones son las siguientes:

|P1 (1.5) − f (1.5)| ≈ 1.53 × 10−3 , |P2 (1.5) − f (1.5)| ≈ 5.42 × 10−4 , | Pˆ 2 (1.5) − f (1.5)| ≈ 6.44 × 10−4 , |P3 (1.5) − f (1.5)| ≈ 2.5 × 10−6 , | Pˆ 3 (1.5) − f (1.5)| ≈ 1.50 × 10−5 , |P4 (1.5) − f (1.5)| ≈ 7.7 × 10−6 . Aunque P3(1.5) es la aproximación más precisa, si no conocemos el valor real de f (1.5), aceptaríamos P4(1.5) como la mejor aproximación ya que incluye la mayor cantidad de datos sobre la función. El término del error de Lagrange derivado del teorema 3.3 no se puede aplicar aquí porque no conocemos la cuarta derivada de f. Por desgracia, este casi siempre es el caso.

Método de Neville 8QDGLÀFXOWDGSUiFWLFDFRQODLQWHUSRODFLyQGH/DJUDQJHHVTXHHOWpUPLQRGHOHUURUHVGLItFLO de aplicar, por lo que el grado del polinomio que se necesita para la precisión deseada en general se desconoce hasta que se realizan los cálculos. Una práctica común es calcular los resultados dados a partir de diferentes polinomios hasta que se obtiene el acuerdo apropiado, como se hizo en la ilustración anterior. Sin embargo, el trabajo efectuado al calcular la aproximación con el segundo polinomio no disminuye el trabajo necesario para calcular la tercera aproximación, ni la cuarta aproximación es fácil de obtener una vez que se conoce la tercera aproximación y así sucesivamente. Ahora, derivaremos estos polinomios de aproximación de una manera que use los cálculos previos para una mayor ventaja. Definición 3.4

Sea fXQDIXQFLyQGHÀQLGDHQx0 , x1 , x2 , . . . , xn y suponga que m 1 , m 2 , . . . , m k son k enteros diferentes, con 0 ≤ m i ≤ n para cada i. El polinomio de Lagrange que concuerda con f(x) en los puntos k xm 1 , xm 2 , . . . , xm k se denota Pm 1 ,m 2 ,... ,m k (x).

88

CAPÍTULO 3

Interpolación y aproximación polinomial

Ejemplo 1

Suponga que x0 5 1, x1 5 2, x2 5 3, x3 5 4, x4 5 6 y f(x) = ex. Determine el polinomio de interpolación que se denota P1,2,4(x) y use este polinomio para aproximar f(5). Solución

Éste es el polinomio de Lagrange que concuerda con f(x) en x1 5 2, x2 5 3 y x4 5 6. Por lo tanto,

P1,2,4 (x) =

(x − 3)(x − 6) 2 (x − 2)(x − 6) 3 (x − 2)(x − 3) 6 e + e + e . (2 − 3)(2 − 6) (3 − 2)(3 − 6) (6 − 2)(6 − 3)

por lo que, f (5) ≈ P(5) =

(5 − 3)(5 − 6) 2 (5 − 2)(5 − 6) 3 (5 − 2)(5 − 3) 6 e + e + e (2 − 3)(2 − 6) (3 − 2)(3 − 6) (6 − 2)(6 − 3)

1 1 = − e2 + e3 + e6 ≈ 218.105. 2 2 El siguiente resultado describe un método para generar de forma recursiva las aproximaciones del polinomio de Lagrange. Teorema 3.5

Sea f GHÀQLGDHQx0, x1, 7, xk y sean xj y xi dos números distintos en este conjunto. Entonces

P(x) =

(x − x j )P0,1,... , j−1, j+1,... ,k (x) − (x − xi )P0,1,... ,i−1,i+1,... ,k (x) (xi − x j )

es el k-ésimo polinomio de Lagrange que interpola f en los puntos k 1 1 x0, x1, 7, xk.

ˆ ≡ P0,1,... , Para la facilidad de la notación, sea Q ≡ P0,1,... ,i−1,i+1,... ,k y Q ˆ son polinomios de grado k 2 1 o menos, P(x) es de j−1, j+1,... ,k . Puesto que Q(x) y Q(x) grado máximo k. ˆ i ) = f (xi ) implica que Primero, observe que Q(x Demostración

P(xi ) =

ˆ i ) − (xi − xi )Q(xi ) (xi − x j ) Q(x (xi − x j ) f (xi ) = f (xi ). = xi − x j (xi − x j )

Similarmente, como Q(x j ) = f (x j ), tenemos que P(x j ) = f (x j ). ˆ r ) = f (xr ). Por lo tanto, Además, si 0 ≤ r ≤ k y r no es i ni j, entonces Q(xr ) = Q(x

P(xr ) =

ˆ r ) − (xr − xi )Q(xr ) (xr − x j ) Q(x (xi − x j ) f (xr ) = f (xr ). = xi − x j (xi − x j )

3HURSRUGHÀQLFLyQP0,1,... ,k (x) es el único polinomio de grado máximo k que concuerda con f en x0 , x1 , . . . , xk . Por lo tanto P ≡ P0,1,... ,k . El teorema 3.5 implica que los polinomios de interpolación pueden generarse de manera recursiva. Por ejemplo, tenemos

P0,1 = P0,1,2 =

1 [(x − x0 )P1 + (x − x1 )P0 ], x1 − x0

P1,2 =

1 [(x − x1 )P2 + (x − x2 )P1 ], x2 − x1

1 [(x − x0 )P1,2 + (x − x2 )P0,1 ], x2 − x0

y así sucesivamente. Estos se generan de la manera que se muestra en la tabla 3.3, donde FDGDÀODVHFRPSOHWDDQWHVGHTXHODVÀODVVXFHVLYDVFRPLHQFHQ

3.2

Tabla 3.3

Eric Harold Neville (1889–1961) DSRUWyHVWDPRGLÀFDFLyQGH la fórmula de Lagrange en un DUWtFXORSXEOLFDGRHQ>1@

x0 x1 x2 x3 x4

P0 P1 P2 P3 P4

P0,1 P1,2 P2,3 P3,4

P0,1,2 P1,2,3 P2,3,4

P0,1,2,3 P1,2,3,4

Aproximación de datos y método de Neville

89

P0,1,2,3,4

El procedimiento que usa el resultado del teorema 3.5 para generar recursivamente las aproximaciones de polinomios de interpolación recibe el nombre de método de Neville. La notación P que se usa en la tabla 3.3 es pesada debido al número de subíndices que se utilizan para representar las entradas. Observe, sin embargo, que mientras se construye un arreglo, sólo se necesitan dos subíndices. El procedimiento hacia abajo en la tabla corresponde al uso consecutivo de los puntos xi con una i más grande, y el procedimiento hacia la derecha corresponde al incremento del grado del polinomio de interpolación. Puesto que los puntos aparecen de manera consecutiva en cada entrada, necesitamos describir sólo un punto de inicio y el número de puntos adicionales que se usan en la construcción de la aproximación. Para evitar los múltiples índices, dejamos que Qi,j (x) para 0 ≤ j ≤ i, denote el polinomio de interpolación de grado j en los números (j + 1) xi− j , xi− j+1 , . . . , xi−1 , xi ; es decir

Q i, j = Pi− j,i− j+1,... ,i−1,i . Usando esta notación obtenemos el arreglo de notación Q en la tabla 3.4.

Tabla 3.4

Ejemplo 2

Tabla 3.5 x

f (x)

1.0 1.3 1.6 1.9 2.2

0.7651977 0.6200860 0.4554022 0.2818186 0.1103623

x0 x1 x2 x3 x4

P0 P1 P2 P3 P4

= = = = =

Q 0,0 Q 1,0 Q 2,0 Q 3,0 Q 4,0

P0,1 P1,2 P2,3 P3,4

= = = =

Q 1,1 Q 2,1 Q 3,1 Q 4,1

P0,1,2 = Q 2,2 P1,2,3 = Q 3,2 P2,3,4 = Q 4,2

P0,1,2,3 = Q 3,3 P1,2,3,4 = Q 4,3

P0,1,2,3,4 = Q 4,4

Los valores de diferentes polinomios de interpolación en x 5 1.5 se obtuvieron en la ilustración al inicio de esta sección usando los datos que se muestran en la tabla 3.5. Aplique el método de Neville a los datos mediante la construcción de una tabla recursiva de la forma que se observa en la tabla 3.4. Solución

Sea x0 5 1.0, x1 5 1.3, x2 5 1.6, x3 5 1.9 y x4 = 2.2, entonces Q0,0 5 f (1.0), Q1,0 5 f(1.3), Q2,0 5 f(1.6), Q3,0 5 f(1.9) y Q4,0 5 f(2.2). Estos son los cinco polinomios de grado cero (constantes) que aproximan f(1.5) y son iguales a los datos que se proporcionan en la tabla 3.5. Al calcular la aproximación de primer grado Q1,1 (1.5) obtenemos

Q 1,1 (1.5) =

(x − x0 )Q 1,0 − (x − x1 )Q 0,0 x1 − x0

(1.5 − 1.0)Q 1,0 − (1.5 − 1.3)Q 0,0 1.3 − 1.0 0.5(0.6200860) − 0.2(0.7651977) = 0.5233449. = 0.3 =

De igual forma, Q 2,1 (1.5) =

(1.5 − 1.3)(0.4554022) − (1.5 − 1.6)(0.6200860) = 0.5102968, 1.6 − 1.3

Q 3,1 (1.5) = 0.5132634,

y

Q 4,1 (1.5) = 0.5104270.

90

CAPÍTULO 3

Interpolación y aproximación polinomial

Se espera que la mejor aproximación lineal sea Q2,1 porque 1.5 se encuentra entre x1 5 1.3 y x2 5 1.6. De manera similar, las aproximaciones usando polinomios de grado superior están dadas por

Q 2,2 (1.5) =

(1.5 − 1.0)(0.5102968) − (1.5 − 1.6)(0.5233449) = 0.5124715, 1.6 − 1.0

Q 3,2 (1.5) = 0.5112857,

Q 4,2 (1.5) = 0.5137361.

y

Las aproximaciones de grado superior se generan de una manera similar y se muestran en la tabla 3.6.

Tabla 3.6

1.0 1.3 1.6 1.9 2.2

0.7651977 0.6200860 0.4554022 0.2818186 0.1103623

0.5233449 0.5102968 0.5132634 0.5104270

0.5124715 0.5112857 0.5137361

0.5118127 0.5118302

0.5118200

Si la última aproximación Q4,4 QRIXHVXÀFLHQWHPHQWHSUHFLVDVHUtDSRVLEOHVHOHFFLRQDU otro nodo x5\DxDGLURWUDÀODDODWDEOD

x5

Q 5,0

Q 5,1

Q 5,2

Q 5,3

Q 5,4

Q 5,5 .

Entonces Q4,4, Q5,4 y Q5,5 podrían compararse para determinar la precisión posterior. /DIXQFLyQHQHOHMHPSORHVODIXQFLyQGH%HVVHOGHSULPHUDFODVHGHRUGHQFHURFX\R valor en 2.5 es 2\ODVLJXLHQWHÀODGHDSUR[LPDFLRQHVSDUDf(1.5) es

2.5

− 0.0483838

0.4807699

0.5301984

0.5119070

0.5118430

0.5118277.

La última nueva entrada, 0.5118277, es correcta para siete lugares decimales. Ejemplo 3

Tabla 3.7 i

xi

ln xi

0 1 2

2.0 2.2 2.3

0.6931 0.7885 0.8329

La tabla 3.7 lista los valores de f(x) 5 ln x precisos para los lugares dados. Use el método de Neville y la aritmética de redondeo de cuatro dígitos para aproximar f (2.1) 5 ln 2.1 al completar la tabla de Neville. Solución Puesto que x 2 x0 5 0.1, x 2 x1 5 20.1 y x 2 x2 5 20.2, tenemos Q0,0 5 0.6931, Q1,0 5 0.7885 y Q2,0 5 0.8329,

Q 1,1 =

1 0.1482 = 0.7410 [(0.1)0.7885 − (−0.1)0.6931] = 0.2 0.2

Q 2,1 =

1 0.07441 = 0.7441. [(−0.1)0.8329 − (−0.2)0.7885] = 0.1 0.1

y

La aproximación final que podemos obtener a partir de estos datos es Q 2,1 =

1 0.2276 = 0.7420. [(0.1)0.7441 − (−0.2)0.7410] = 0.3 0.3

Estos valores se muestran en la tabla 3.8.

Tabla 3.8

i

xi

x − xi

Q i0

Q i1

Q i2

0 1 2

2.0 2.2 2.3

0.1 −0.1 −0.2

0.6931 0.7885 0.8329

0.7410 0.7441

0.7420

3.3

Diferencias divididas

91

En el ejemplo anterior, tenemos f(2.1) 5 ln 2.1 5 0.7419 para cuatro lugares decimales, por lo que el error absoluto es

| f (2.1) − P2 (2.1)| = |0.7419 − 0.7420| = 10−4 . Sin embargo, f (x) = 1/x, f (x) = −1/x 2 , y f (x) = 2/x 3, por lo que la fórmula de error de Lagrange (3.3) en el teorema 3.3 nos da la cota del de error

| f (2.1) − P2 (2.1)| = =

f (ξ(2.1)) (x − x 0 )(x − x1 )(x − x2 ) 3! 1 0.002 (0.1)(−0.1)(−0.2) ≤ = 8.3 × 10−5 . 3 3(2)3 3 (ξ(2.1))

Observe que el error real, 1024, excede la cota del error, 8.3 × 10−5. Esta aparente conWUDGLFFLyQHVXQDFRQVHFXHQFLDGHORVFiOFXORVGHGtJLWRVÀQLWRV1RVRWURVXVDPRVODDULWmética de redondeo de cuatro dígitos, y la fórmula del error de Lagrange (3.3) supone la DULWPpWLFDGHGtJLWRVLQÀQLWRV(VWRFDXVyTXHQXHVWURVHUURUHVUHDOHVH[FHGLHUDQHOFiOFXOR de error teórico. • Recuerde: No puede esperar mayor precisión de la proporcionada por la aritmética. (ODOJRULWPRFRQVWUX\HSRUÀODVODVHQWUDGDVHQHOPpWRGRGH1HYLOOH ALGORITMO

3.1

Interpolación iterada de Neville Para evaluar el polinomio de interpolación P en los diferentes números n 1 1, x0, 7, xn en el número x para la función f :

ENTRADA números x, x0 , x1 , . . . , xn ; valores f (x0 ), f (x1 ), . . . , f (xn ) como la primera columna Q 0,0 , Q 1,0 , . . . , Q n,0 de Q. SALIDA

la tabla Q con P(x) = Q n,n .

Paso 1 Para i = 1, 2, . . . , n para j = 1, 2, . . . , i haga Q i, j =

(x − xi− j )Q i, j−1 − (x − xi )Q i−1, j−1 . xi − xi− j

Paso 2 SALIDA (Q); PARE. La sección Conjunto de ejercicios 3.2 está disponible en línea. Encuentre la ruta de en las páginas preliminares.

3.3 Diferencias divididas La interpolación iterada se usó en la sección previa para generar sucesivamente aproximaFLRQHVSROLQRPLDOHVGHJUDGRVXSHULRUHQXQSXQWRHVSHFtÀFR/RVPpWRGRVGHGLIHUHQFLD dividida que se presentan en esta sección se usan para generar sucesivamente los polinomios en sí mismos.

Diferencias divididas Suponga que Pn(x) es el enésimo polinomio de interpolación que concuerda con la función f en los diferentes números x0, x1, 7, xn. A pesar de que este polinomio es único, existen re-

92

CAPÍTULO 3

Interpolación y aproximación polinomial

presentaciones algebraicas que son útiles en ciertas situaciones. Las diferencias divididas de f respecto a x0, x1, 7, xn se usan para expresar Pn(x) en la forma

Pn (x) = a0 + a1 (x − x0 ) + a2 (x − x0 )(x − x1 ) + · · · + an (x − x0 ) · · · (x − xn−1 ), (3.5) para constantes apropiadas a0, a1, 7, an. Para determinar la primera de estas constantes, a0, observe que si Pn(x) se escribe en la forma de la ecuación (3.5), entonces evaluando Pn(x) en x0 queda sólo el término constante a0; es decir,

a0 = Pn (x0 ) = f (x0 ). Como en muchas áreas, Isaac Newton es prominente en el estudio de ecuaciones de diferencia. Desarrolló fórmulas de interpolación desde 1675, usando su notación en tablas de diferencias. Adoptó un enfoque muy general hacia las fórmulas de diferencias, por lo que los ejemplos explícitos que produjo, incluyendo las fórmulas de Lagrange, a menudo son conocidas con otros nombres.

Similarmente, cuando P(x) se evalúa en x1, los únicos términos diferentes de cero en la evaluación de Pn(x1) son los términos constante y lineal,

f (x0 ) + a1 (x1 − x0 ) = Pn (x1 ) = f (x1 ); por lo que a1 =

f (x1 ) − f (x0 ) . x1 − x0

(3.6)

Ahora presentaremos la notación de diferencias divididas, que se relaciona con la notación 2 de Aitkens que se usó en la sección 2.5. La ceroésima diferencia dividida de la función f respecto a xi, denotada f >xi], es simplemente el valor de f en xi:

f [xi ] = f (xi ).

(3.7)

/DVGLIHUHQFLDVGLYLGLGDVUHVWDQWHVVHGHÀQHQGHPDQHUDUHFXUVLYDODprimera diferencia dividida de f respecto a xi y xi+1 se denota f [xi , xi+1 ] \VHGHÀQHFRPR

f [xi , xi+1 ] =

f [xi+1 ] − f [xi ] . xi+1 − xi

(3.8)

La segunda diferencia dividida, f [xi , xi+1 , xi+2 ], VHGHÀQHFRPR

f [xi , xi+1 , xi+2 ] =

f [xi+1 , xi+2 ] − f [xi , xi+1 ] . xi+2 − xi

De igual forma, después de que las (k 21) -ésimas diferencias divididas,

f [xi , xi+1 , xi+2 , . . . , xi+k−1 ] y

f [xi+1 , xi+2 , . . . , xi+k−1 , xi+k ],

se han determinado, la k-ésima diferencia dividida relativa a xi , xi+1 , xi+2 , . . . , xi+k es

f [xi , xi+1 , . . . , xi+k−1 , xi+k ] =

f [xi+1 , xi+2 , . . . , xi+k ] − f [xi , xi+1 , . . . , xi+k−1 ] . xi+k − xi

(3.9)

El proceso termina con la única enésima diferencia dividida,

f [x0 , x1 , . . . , xn ] =

f [x1 , x2 , . . . , xn ] − f [x0 , x1 , . . . , xn−1 ] . xn − x0

Debido a la ecuación (3.6), podemos escribir a1 = f [x0 , x1 ], justo cuando a0 se puede expresar como a0 = f (x0 ) = f [x0 ]. Por lo tanto, el polinomio de interpolación en la ecuación (3.5) es

Pn (x) = f [x0 ] + f [x0 , x1 ](x − x0 ) + a2 (x − x0 )(x − x1 ) + · · · + an (x − x0 )(x − x1 ) · · · (x − xn−1 ).

Análisis numérico, 10a. ed., se escribió para que los estudiantes de ingeniería, matemáticas, ciencias de la computación puedan usarlo en los cursos sobre la teoría y la aplicación de técnicas de aproximación numérica. Prácticamente todos los conceptos en el texto están ilustrados con un ejemplo y contiene más de 2 500 ejercicios probados en clase que van desde aplicaciones fundamentales de métodos y algoritmos hasta generalizaciones y extensiones de la teoría. Además, los conjuntos de ejercicios incluyen varios problemas aplicados de diversas áreas de la ingeniería, así como de la física, la informática, la biología y las ciencias económicas y sociales. Las aplicaciones, seleccionadas de forma clara y concisa, demuestran la manera en la que las técnicas numéricas se aplican en situaciones de la vida real.

Visite nuestro sitio en http://latinoamerica.cengage.com