Capítulo)3.)Derivación)e)integración)numérica)aproximada)
de)funciones)
Métodos)Numéricos)
Carlos)Beltrán)Álvarez)
Departamento*de*Matemá.cas,*Estadís.ca*y*
Computación*
Este*tema*se*publica*bajo*Licencia:*
Crea.ve*Commons*BYANCASA*4.0*
Prepara tus exámenes
Prepara tus exámenes y mejora tus resultados gracias a la gran cantidad de recursos disponibles en Docsity
Consigue puntos base para descargar
Gana puntos ayudando a otros estudiantes o consíguelos activando un Plan Premium
Orientación Universidad
Prepara tus exámenes
Prepara tus exámenes y mejora tus resultados gracias a la gran cantidad de recursos disponibles en Docsity
Busca documentos
Prepara tus exámenes con los documentos que comparten otros estudiantes como tú en Docsity
Busca tu universidad
Encuentra los documentos específicos para los exámenes de tu universidad
Video Cursos
Estudia con lecciones y exámenes resueltos basados en los programas académicos de las mejores universidades
Quiz
Responde a preguntas de exámenes reales y pon a prueba tu preparación
Consigue puntos base para descargar
Gana puntos ayudando a otros estudiantes o consíguelos activando un Plan Premium
Comunidad
Pregúntale a la comunidad
Pide ayuda a la comunidad y resuelve tus dudas de estudio
Ebooks gratuitos
¡Nuestros e-books salva-estudiantes!
Descarga nuestras guías gratuitas sobre técnicas de estudio, métodos para controlar la ansiedad y consejos para la tesis preparadas por los tutores de Docsity
Método de la derivada de tres puntos y aproximación de derivadas y integrales, Apuntes de Informática
El método de la derivada de tres puntos y su aplicación para aproximar derivadas y integrales. Se incluyen formulas para derivadas y integrales de diferentes funciones, como polinomios, exponenciales, trigonométricas y combinadas. Se utiliza el teorema de Taylor para decir más sobre las aproximaciones. Además, se presenta una demostración independiente del teorema de Simpson.
Tipo: Apuntes
2019/2020
1 / 31
Esta página no es visible en la vista previa
¡No te pierdas las partes importantes!
Documentos relacionados
Vista previa parcial del texto
¡Descarga Método de la derivada de tres puntos y aproximación de derivadas y integrales y más Apuntes en PDF de Informática solo en Docsity!
Capítulo 3. Derivación e integración numérica aproximada
de funciones
Métodos Numéricos
Carlos Beltrán Álvarez
Departamento de Matemá.cas, Estadís.ca y
Computación
Este tema se publica bajo Licencia:
Crea.ve Commons BY-‐NC-‐SA 4.
Cap´ıtulo 3
Derivaci´on e integraci´on num´erica
aproximada de funciones
En este cap´ıtulo nos enfrentamos a uno de los problemas fundamentales del c´alculo (encontrar derivadas e integrales de funciones dadas) desde la nueva perspectiva de los m´etodos num´ericos. Es habitual que en los cursos iniciales de c´alculo uno aprenda a derivar b´asicamente cualquier funci´on escrita mediante una f´ormula, y en efecto utilizando las reglas de la suma, el producto, el cociente, la regla de la cadena y la regla de la derivada de la funci´on inversa se obtienen pr´acticamente todas las derivadas que uno desee de una funci´on dada. Sin embargo, hay situaciones en las que este procedimiento, tan habitual para nosotros, no es factible. Algunas de estas situaciones son:
i) La funci´on que deseamos derivar tiene una f´ormula extremadamente complicada, por ejemplo deseamos calcular la segunda derivada en x = 0 de una funci´on que contiene varios productos y composiciones de funciones elementales y resulta muy trabajoso hacer todo el c´alculo simb´olico para luego hacer x = 0 y obtener un simple n´umero. Es natural preguntarse si no habr´a un m´etodo menos trabajoso teniendo en cuenta que la informaci´on requerida es tambi´en menor (solo queremos la derivada en un punto, no una f´ormula general).
ii) La funci´on que deseamos derivar se calcula mediante un programa complejo, por ejemplo un programa que en input t devuelve el valor de un cierto trazador c´ubico en el punto t. En este caso no tenemos una f´ormula, lo ´unico que podemos es generar, usando el programa, los valores de la funci´on en distintos puntos. Otra situaci´on m´as realista viene dada por el ejemplo III en la introducci´on del Cap´ıtulo 4
iii) La funci´on existe pero, como ya nos hemos encontrado en el cap´ıtulo de interpolaci´on, solo podemos conocerla en algunos puntos concretos con medidas experimentales. Por ejemplo, observamos en dos d´ıas consecutivos el planeta Marte a la misma hora de la noche, de forma que las estrellas fijas pr´acticamente no se han movido pero el planeta s´ı. A partir de estas observaciones, ¿c´omo deducimos el valor de la velocidad angular del planeta? Solo tenemos dos datos, pero veremos que son suficientes si le hacemos algunas concesiones a la precisi´on.
El caso de los problemas de integraci´on es m´as familiar para todos: ¿qui´en no se ha encontrado con la necesidad –acad´emica o cient´ıfica– de calcular una integral, siendo imposible resolverla mediante las t´ecnicas habituales? En el coraz´on de este problema est´a el conocido resultado (escrito aqu´ı en su forma m´as elemental):
3.1. C ALCULO APROXIMADO DE DERIVADAS´ 41
Teorema 3.1.1 (Derivada hacia adelante) Sea f ∈ C^1 [a, b] tal que existe f ′′^ en (a, b). Sean x 0 , h ∈ R, h 6 = 0, tales que x 0 , x 0 + h ∈ [a, b]. Entonces,
f ′(x 0 ) =
f (x 0 + h) − f (x 0 ) h
h 2
f ′′(ζ),
para alg´un ζ ∈ (a, b) entre x 0 y x 0 + h (que depende de x 0 y de h). N´otese que puede ser h < 0.
Demostraci´on. Por el Teorema 1.3.1 tenemos f (x 0 + h) = f (x 0 ) + f ′(x 0 )h + h^2 f ′′(ζ)/2 para alg´un ζ entre x 0 y x 0 + h. Por lo tanto,
f ′(x 0 ) −
f (x 0 + h) − f (x 0 ) h
h 2
f ′′(ζ),
como quer´ıamos. N´otese que el Teorema 3.1.1 solo nos permite acotar el error si la derivada la calculamos en un nodo del polinomio interpolante. Nada nos impide utilizar la f´ormula en otro caso, si estamos dispuestos a renunciar a tener una cota para el error.
Ejemplo 3.1.2 Consideremos la funci´on f (x) = e(log^ x) 2 definida para x > 0. Calculemos la derivada de f en x = 1 utilizando el m´etodo de derivada hacia delante con varios valores de h. Con h = 0, 1 , tenemos
f ′(1) ≈ f (1,1) − f (1) 0 , 1
y de acuerdo con el Teorema 3.1.1 el error cometido es de la forma − 0 , 05 f ′′(ζ) para alg´un ζ ∈ (1, 1 ,1). En situaciones reales no siempre podemos acotar la segunda derivada, pero en este caso lo hacemos por motivos de ilustraci´on. Tenemos
f ′(x) = 2
log x x
e(log^ x)
2 , f ′′(x) =
log x x
1 − log x x^2
e(log^ x)
2 .
Vemos que f ′′^ es positiva en [1, 1 ,1] y podemos acotarla f´acilmente en este intervalo: para x ∈ [1, 1 ,1] tenemos
0 ≤ f (x) ≤
log 1, 1 1 , 1
e(log 1,1)
2 ≤ 2 , 0789.
Con ello, podemos asegurar que nuestro c´alculo aproximado de f ′(1) comete un error en el intervalo [− 0 , 05 · 2 , 0789 , 0], esto es el error cometido est´a entre − 0 , 1039 y 0. Como en este ejemplo podemos calcular el valor exacto de la derivada f ′(1) = 0, comprobamos que el error cometido (igual a − 0 , 0913 ) est´a en efecto en el intervalo de error que nos da el Teorema 3.1.1.
Derivada de tres puntos y derivada central (o del punto medio)
El siguiente paso en nuestro m´etodo general consiste en tomar el polinomio interpolante de grado 2 en 3 nodos, y luego calcular la derivada de la funci´on en uno de esos nodos. Tomemos pues x 0 , x 1 , x 2 ∈ [a, b] en los que el valor de la funci´on es conocido. El polinomio interpolante del Teorema 2.1.2 es entonces
p(x) = f (x 0 ) (x − x 1 )(x − x 2 ) (x 0 − x 1 )(x 0 − x 2 )
f (x 1 ) (x − x 0 )(x − x 2 ) (x 1 − x 0 )(x 1 − x 2 )
f (x 2 ) (x − x 0 )(x − x 1 ) (x 2 − x 0 )(x 2 − x 1 )
42 CAP´ITULO 3. DERIVACI ON E INTEGRACI ´ ON NUM ´ ERICA APROXIMADA DE FUNCIONES´
Por lo tanto, tenemos:
f ′(x) ≈ p′(x) = f (x 0 )
2 x − x 1 − x 2 (x 0 − x 1 )(x 0 − x 2 )
- f (x 1 )
2 x − x 0 − x 2 (x 1 − x 0 )(x 1 − x 2 )
- f (x 2 )
2 x − x 0 − x 1 (x 2 − x 0 )(x 2 − x 1 )
que es la llamada f´ormula de tres puntos. Poniendo ahora x 1 = x 0 + h 1 , x 2 = x 0 + h 2 tenemos:
f ′(x 0 ) ≈ −f (x 0 )
h 1 + h 2 h 1 h 2
− f (x 0 + h 1 )
h 2 h 1 (h 1 − h 2 )
− f (x 0 + h 2 )
h 1 h 2 (h 2 − h 1 )
De nuevo mostramos una expresi´on para el error:
Teorema 3.1.3 Sea f ∈ C^3 [a, b]. Sean x 0 , h 1 , h 2 ∈ R, h 1 , h 2 6 = 0, h 1 6 = h 2 , tales que x 0 , x 0 + h 1 , x 0 + h 2 ∈ [a, b]. Entonces,
f ′(x 0 ) = −f (x 0 )
h 1 + h 2 h 1 h 2
− f (x 0 + h 1 )
h 2 h 1 (h 1 − h 2 )
− f (x 0 + h 2 )
h 1 h 2 (h 2 − h 1 )
h 1 h 2 6
f ′′′(ζ),
para alg´un ζ ∈ (a, b) entre x 0 , x 0 + h 1 y x 0 + h 2 (que depende de x 0 , h 1 , h 2 ). N´otese que puede ser h 1 < 0 y/o h 2 < 0.
Demostraci´on. Por el Teorema 2.1.6, en las condiciones del teorema tenemos:
f (x) = p(x) +
f ′′′(ζx) 6
(x − x 0 )(x − x 1 )(x − x 2 ),
donde para cada x ∈ [a, b] el punto ζx est´a en (a, b). Despejando el t´ermino con la tercera derivada tenemos:
f ′′′(ζx) =
6(f (x) − p(x)) (x − x 0 )(x − x 1 )(x − x 2 )
6(f (x) − f (x 0 )) (x − x 0 )(x − x 1 )(x − x 2 )
6(p(x 0 ) − p(x)) (x − x 0 )(x − x 1 )(x − x 2 )
que tiene l´ımite cuando x → x 0 igual a
6 (x − x 1 )(x − x 2 )
(f ′(x 0 ) − p′(x 0 )). (3.2)
Por compacidad del intervalo [a, b], existe un punto de acumulaci´on ζ del conjunto {ζx : x ∈ [a, b]} tal que f ′′′(ζ) alcanza ex´actamente el valor (3.2). Despejando f ′(x 0 ), obtenemos:
f ′(x 0 ) = p′(x 0 ) + f ′′′(ζ) 6
(x 0 − x 1 )(x 0 − x 2 ).
El resultado se sigue notando que x 0 − x 1 = −h 1 , x 0 − x 2 = −h 2. La f´ormula de la derivada de tres puntos adquiere una forma particularmente sencilla en tres casos:
Si h 2 = h, h 1 = −h entonces el Teorema 3.1.3 dice:
f ′(x 0 ) =
f (x 0 + h) − f (x 0 − h) 2 h
h^2 6
f ′′′(ζ), para alg´un ζ ∈ (x 0 − h, x 0 + h), (3.3)
que se conoce como f´ormula de la derivada central o del punto medio, y es tiene una buena precisi´on con solo dos evaluaciones de la funci´on.
44 CAP´ITULO 3. DERIVACI ON E INTEGRACI ´ ON NUM ´ ERICA APROXIMADA DE FUNCIONES´
Derivada de cinco puntos
De igual manera que hemos hecho en la secci´on anterior con un polinomio cuadr´atico interpolado en 3 puntos, podemos hacerlo para orden superior. Obviamente los c´alculos se complican, pero de entre todas las f´ormulas que se pueden obtener hay una particularmente sencilla que escribimos sin demostraci´on. Se obtiene interpolando la funci´on en los 5 puntos dados por x 0 , x 0 ± h, x 0 ± 2 h y aproximando la derivada de la funci´on por la derivada del polinomio interpolante. La f´ormula obtenida, que presentamos a continuaci´on, se conoce como f´ormula de la derivada de cinco puntos, y es en general m´as precisa que la de la derivada central, aunque exige evaluar la funci´on en 4 puntos en lugar de solo en 2.
Teorema 3.1.6 Sea f ∈ C^5 [a, b], y sean x 0 ∈ [a, b], h > 0 tales que x 0 ± h, x 0 ± 2 h ∈ [a, b]. Entonces,
f ′(x 0 ) =
f (x 0 − 2 h) − 8 f (x 0 − h) + 8f (x 0 + h) − f (x 0 + 2h) 12 h
h^4 30
f (5)(ζ),
para alg´un ζ ∈ (x 0 − 2 h, x 0 + 2h).
Ejemplo 3.1.7 Calculemos de forma aproximada la derivada de la funci´on f (x) = e(log^ x)
2 del Ejemplo 3.1.2, utilizando la f´ormula de la derivada de 5 puntos con h = 0, 1 :
f ′(1) ≈
f (0,8) − 8 f (0,9) + 8f (1,1) − f (1,2) 1 , 2
que es una aproximaci´on de f ′(1) = 0 mejor que la proporcionada por la derivada hacia delante y mejor tambi´en que la proporcionada por la f´ormula de la derivada central.
3.1.2. F´ormulas para la segunda derivada
Al igual que en el caso de la primera derivada, para calcular la segunda derivada tomamos el polinomio interpo- lante de la funci´on en una colecci´on de nodos y aproximamos la segunda derivada de la funci´on por la segunda derivada del polinomio. Tomando un polinomio cuadr´atico como en (3.1) tenemos la aproximaci´on:
f ′′(x 0 ) ≈ p′′(x 0 ) =
2 f (x 0 ) (x 0 − x 1 )(x 0 − x 2 )
2 f (x 1 ) (x 1 − x 0 )(x 1 − x 2 )
2 f (x 2 ) (x 2 − x 0 )(x 2 − x 1 )
donde debemos se˜nalar que la aproximaci´on de la segunda derivada no depende del punto, esto es, si elegimos por ejemplo x 0 = 1, x 1 = 1, 1 , x 2 = 0,9 concretos y usamos la f´ormula para aproximar f ′′(1) obtenemos el mismo n´umero que si la usamos para aproximar por ejemplo f ′′(8), lo que suena un poco absurdo. Sin embargo, no es absurdo: la precisi´on de nuestro c´alculo s´ı que depender´a del punto. Llamando como antes x 1 = x 0 + h 1 , x 2 = x 0 + h 2 concluimos:
f ′′(x 0 ) ≈
2 f (x 0 ) h 1 h 2
2 f (x 0 + h 1 ) h 1 (h 1 − h 2 )
2 f (x 0 + h 2 ) h 2 (h 2 − h 1 )
De forma similar al caso de la primera derivada, tenemos una f´ormula m´as sencilla suponiendo que x 1 = x 0 − h y x 2 = x 0 + h
Teorema 3.1.8 Sea f ∈ C^4 [a, b]. Sean x 0 , h ∈ R, h > 0 , tales que x 0 ± h ∈ [a, b]. Entonces,
f ′′(x 0 ) =
f (x 0 + h) + f (x 0 − h) − 2 f (x 0 ) h^2
h^2 12
f (4)(ζ),
para alg´un ζ ∈ (x 0 − h, x 0 + h).
3.1. C ALCULO APROXIMADO DE DERIVADAS´ 45
Demostraci´on. Por el Teorema 1.3.1 tenemos
f (x 0 + h) + f (x 0 − h) − 2 f (x 0 ) h^2
h^2 f ′′(x 0 ) + h
4 24 (f^
(4)(ζ 1 ) + f (4)(ζ 2 )) h^2
para algunas ζ 1 , ζ 2 ∈ (x 0 − h, x 0 + h). El resultado se sigue usando el Teorema 1.1.7. Para aproximar la segunda derivada de f : [a, b] → R en los extremos del intervalo, seguimos el mismo procedi- miento obteniendo las f´ormulas
f ′′(a) ≈ f (a + 2h) + f (a) − 2 f (a + h) h^2
, f ′′(a) ≈ f (a − 2 h) + f (a) − 2 f (a − h) h^2
Ejemplo 3.1.9 Utilizemos la f´ormula del Teorema 3.1.8 para calcular la segunda derivada de f (x) = e(log^ x)
2
en x = 1 con h = 0, 1 :
f ′′(x) ≈
f (1,1) + f (0,9) − 2 f (1) 0 , 01
mientras que el valor exacto es f ′′(1) = 2.
Ejemplo 3.1.10 Consideremos de nuevo los datos de la Figura 2.3, y calculemos de forma aproximada la segunda derivada de la presi´on en funci´on de la temperatura en T = 2, 7. Para ello usaremos la f´ormula (3.6) con h 1 = − 0 , 4 , h 2 = 0, 2 , obteniendo:
d^2 P dT 2
2 P (2,7)
2 P (2,3)
2 P (2,9)
Vamos a utilizar esta informaci´on de la siguiente manera: en el Ejemplo 3.1.5 hemos calculado de manera aproximada la primera derivada P ′(2,7) ≈ 22 , 8779. Podemos por lo tanto ya aproximar el polinomio de Taylor de segundo orden de P cerca de T = 2, 7 :
q(T ) = P (2,7) + P ′(2,7)(x − 2 .t) +
P ′′(2,7)
(x − 2 ,7)^2 ≈ 13 ,6218 + 22,8779(x − 2 ,7) + 11,9655(x − 2 ,7)^2.
En la Figura 3.1 dibujamos dicho polinomio junto con la tabla de datos.
3.1.3. Adivinando el futuro sin saber casi nada
Con estos resultados para c´alculo aproximado de derivadas y segundas derivadas, y armados del Teorema de Taylor podemos calcular aproximadamente el valor en el futuro inmediato de una funci´on y : [− 2 h, h] → R a partir de su valor actual y(0), y su valor en dos posiciones anteriores equiespaciadas y(−h), y(− 2 h):
y(h) ≈ y(0) + hy′(0) +
h^2 2
y′′(0) ≈ y(0) + h
− 4 y(−h) + 3y(0) + y(− 2 h) 2 h
h^2 2
y(− 2 h) + y(0) − 2 y(−h) h^2
y simplificando: y(h) ≈ y(− 2 h) + 3(y(0) − y(−h)). (3.7)
Dado que la derivada de una funci´on vectorial se calcula componente a componente, esta misma aproximaci´on vale si y : [− 2 h, h] → Rn. Esta observaci´on resulta particularmente interesante si pensamos que la funci´on define la trayectoria de un objeto, r : [− 2 h, h] → R^3. Supongamos que conocemos r(− 2 h), r(−h) y r(0), esto es la
3.1. C ALCULO APROXIMADO DE DERIVADAS´ 47
Se˜nalamos que para utilizar este m´etodo de predicci´on no necesitamos conocer nada de los procesos que subyacen a la funci´on de estudio, esto es, no necesitamos conocer las leyes f´ısicas que rigen los procesos para predecirlos a corto plazo. Los ni˜nos peque˜nos, que por lo general no saben de F´ısica, pueden predecir f´acilmente d´onde caer´a una pelota que se les lanza, con suficiente precisi´on como para alcanzarla en muchas ocasiones. A corto plazo, no hace falta demasiada informaci´on, solo hace falta saber utilizarla! Veamos otro ejemplo de uso.
Ejemplo 3.1.12 Hemos medido la altitud (altura sobre el horizonte) de Saturno en el cielo nocturno tres d´ıas consecutivos a las 23 : 00, a saber, los d´ıas 2 , 3 y 4 de julio de 2017, que fueron tres d´ıas particularmente buenos para la observaci´on astron´omica en Cantabria. Los valores respectivos de la latitud fueron 19 o 05 ′, 19 o 27 ′, 19 o 49 ′. En grados: 19 , 08 o, 19 , 45 o, 19 , 82 o.
Utilizando la f´ormula (3.7) podemos dar una aproximaci´on de la altitud de Saturno a la misma hora del d´ıa 5 de julio (d´ıa en el que este ejemplo se ha escrito). La predicci´on es:
19 ,08 + 3 · (19, 82 − 19 ,45) = 20, 19 o^ = 20o 11 ′.
La altitud real de Saturno a las 23:00 horas de la noche del 5 de julio fue de 20 o 10 , con lo que nuestra predicci´on se equivoc´o por un minuto de arco. N´otese que no hemos necesitado utilizar ning´un conocimiento sobre las caracter´ısticas orbitales de Saturno ni ninguna otra informaci´on para realizar nuestra predicci´on.
3.1.4. C´alculo de derivadas en presencia de errores de medici´on
Mirando a las f´ormulas de los teoremas anteriores para el c´alculo de derivadas, da la impresi´on de que lo que conviene es tomar un valor de h > 0 extremadamente peque˜no, por ejemplo h = 10−^16 o incluso h = 10−^1000. Sin embargo, nuestra intuici´on nos dice que algo fallar´a si lo hacemos as´ı: ser´a muy dif´ıcil distinguir los valores de f (x 0 − h), f (x 0 ) y f (x 0 − h) con lo que seguramente nuestro c´alculo no es muy preciso. De hecho, para calcular la derivada estaremos llevando a cabo dos de las operaciones poco recomendables seg´un la Secci´on 1.4: restar dos n´umeros muy parecidos y dividir por un n´umero muy peque˜no. En esta secci´on le damos un contexto te´orico a esta observaci´on, partiendo de la suposici´on de que las medidas que podemos hacer de la funci´on f tienen un peque˜no error .
Teorema 3.1.13 Supongamos que podemos calcular valores de f ∈ C^3 [a, b] con precisi´on . Entonces, utilizando la f´ormula de la primera derivada central, en el mejor de los casos podemos calcular valores de f ′^ con precisi´on del orden de ^2 /^3.
Demostraci´on. Nuestro c´alculo de la derivada en x 0 ser´a:
f (x 0 + h) + 1 − f (x 0 − h) − 2 2 h
= f ′(x 0 ) + h^2 6
f ′′′(ζ) +
2 h
donde ζ ∈ (x 0 − h, x 0 + h) y lo ´unico que sabemos de 1 , 2 es que son reales de m´odulo a lo m´as . El error en valor absoluto puede por lo tanto alcanzar (dado que los signos de los pueden ser distintos) un valor de
h
h^2 6
|f ′′′(ζ)| ≤
h
h^2 6
M,
donde M es alguna cota para la tercera derivada en [a, b]. Si elegimos el mejor h posible, esto es, el que minimiza el error (derivando e igualando a 0 la expresi´on de arriba), tenemos:
h =
M
48 CAP´ITULO 3. DERIVACI ON E INTEGRACI ´ ON NUM ´ ERICA APROXIMADA DE FUNCIONES´
y por lo tanto la cota que podemos dar para el error es: ( 9 M 8
^2 /^3.
En la demostraci´on anterior no conocemos el valor de M , lo que complica las cosas en la pr´actica pues no podemos elegir el mejor valor de h, pero lo que dice el teorema es que sea cual sea el valor de M (salvo que sea 0 o extremadamente cercano a 0, esto es, salvo que la funci´on que estamos derivando sea pr´acticamente igual a una par´abola), el error que cometemos al calcular la derivada es mucho mayor que el que cometemos al calcular la funci´on. Por ejemplo si calculamos la funci´on con precisi´on de una millon´esima es = 10−^6 pero la derivada la calcularemos ´unicamente con precisi´on del orden de ^2 /^3 = 10−^4 , esto es del orden de una diezmil´esima. En la pr´actica tratamos de elegir valores de h que no sean ni demasiado grandes ni demasiado peque˜nos, y solo una cierta familiaridad con el problema nos permite elegirlos de esa manera. N´otese que dependiendo del problema que estemos tratando un mismo n´umero puede considerarse “grande” o “peque˜no” (por ejemplo, 10 segundos en el movimiento de la Tierra al rededor del Sol es un tiempo extremadamente peque˜no, pero 10 segundos para un objeto en ca´ıda libre en la superficie de la Luna es un tiempo bastante grande). Usando la f´ormula de cinco puntos la situaci´on mejora pero solo un poco. Para ´ordenes superiores la situaci´on es l´ogicamente peor. Queda como tarea para el alumno el calcular el orden de precisi´on de las derivadas obtenidas mediante estos m´etodos, esto es, resultados an´alogos al Teorema 3.1.13 para la f´ormula de cinco puntos, para la f´ormula de la segunda derivada sim´etrica, etc.
Ejemplo 3.1.14 Para ilustrar c´omo el uso de un h dado influye en la aproximaci´on de la derivada, calculamos la derivada de f (x) = ex^ en x = 0 (cuyo valor exacto es f ′(0) = 1) de forma aproximada con distintos valores de h, mostrando el resultado en la Figura 3.2.
3.1.5. Derivadas de orden superior
Si en alg´un momento necesitamos calcular por ejemplo la tercera o la cuarta derivada de una funci´on, podemos utilizar el siguiente m´etodo:
f ′′′(x) ≈
f ′′(x + h) − f ′′(x − h) 2 h
y a continuaci´on la f´ormula del Teorema 3.1.8 para aproximar f ′′(x ± h). Obviamente tenemos m´as alternativas, por ejemplo utilizar
f ′′′(x) ≈
f ′(x + h) + f ′(x − h) − 2 f ′(x) h^2
y entonces usar la f´ormula del Teorema 3.1.3 para aproximar f ′(x ± h). No damos detalles sobre las f´ormulas de error, pero s´ı que indicamos que tal como se ha explicado en la Secci´on 3.1.4 el c´alculo ser´a menos preciso para ´ordenes de derivaci´on superior, y solo lo acometeremos cuando no tengamos otro remedio.
3.1.6. C´alculo de matrices Jacobianas, gradientes, divergencias y Laplacianos
Es frecuente tambi´en tener que calcular matrices Jacobianas y gradientes de funciones. Recordemos que si f : Rn^ → Rm^ entonces la matriz Jacobiana de f en v ∈ Rm^ se define por
Jf (v) =
∂f 1 ∂x 1 · · ·^
∂f 1 ∂xn .. .
∂fm ∂x 1 · · ·^
∂fm ∂xn
∂f ∂x 1 · · ·^
∂f ∂xn
50 CAP´ITULO 3. DERIVACI ON E INTEGRACI ´ ON NUM ´ ERICA APROXIMADA DE FUNCIONES´
for k=1:n % k-esima columna de J: e=zeros(n,1); e(k)=1; % k-esimo vector de la base canonica. J(:,k)=(f(v+he)-f(v-he))/2/h; end
Si deseamos calcular el gradiente de f : Rn^ → R ya sabemos que basta tomar Jf (x)T^. En cuanto a la divergencia de una funci´on vectorial f : Rn^ → Rn, definida por
divf (v) =
∂f 1 ∂x 1
∂fn ∂xn
podemos utilizar un procedimiento similar (calculando cada derivada parcial con el m´etodo dado en (3.8) o directamente tomando la traza de la matriz Jacobiana, un procedimiento m´as costoso en tiempo de ejecuci´on pero el que usamos a continuaci´on):
function div=divergencia(f,v,h) % Calcula la divergencia de f en v, usando la formula central. % Se supone que f va de Rˆn a Rˆn, todos los vectores son columnas. % h es el parametro para calcular la derivada J=matrizjacobiana(f,v,h); div=trace(J);
Ejemplo 3.1.15 Calcula de forma aproximada el gradiente de la funci´on f (x, y, z) = xy cos(z)+exz^ en el punto (x, y, z) = (1, 1 , π), utilizando la f´ormula de la derivada central con paso h = 10−^1. Escribe una estimaci´on en error absoluto del error que cometes en cada coordenada del gradiente. Respuesta: tenemos ∇f (1, 1 , π) ≈ (a, b, c) con:
a = f (1, 1 , 1 , π) − f (0, 9 , 1 , π) 0 , 2
b = f (1, 1 , 1 , π) − f (1, 0 , 9 , 1 , π) 0 , 2
c =
f (1, 1 , 3 ,241592) − f (1, 1 , 3 ,041592) 0 , 2
Estimamos los errores: dado que a = g′(0) con g(t) = f (1 + t, 1 , π), la f´ormula para el valor absoluto del error cometido es (para alg´un ζ ∈ (− 0 , 1 , 0 ,1)):
h^2 6
|g′′′(ζ)| =
π^3 e(1+ζ)π 600
De forma similar, para b y c los errores son respectivamente a lo sumo:
0 600
sin 0,1 + eπ+0,^1 600
Finalmente, damos tambi´en un m´etodo para calcular el Laplaciano de f : Rn^ → R:
∆f (v) =
∂^2 f ∂x^21
∂^2 f ∂x^2 n
3.1. C ALCULO APROXIMADO DE DERIVADAS´ 51
aproximando luego cada derivada segunda por la f´ormula
∂^2 f ∂x^2 j
f (v + hej ) + f (v − hej ) − 2 f (v) h^2
que consiste en aplicar el m´etodo del Teorema 3.1.8. Con ello, un programa de Matlab para calcular el Laplaciano es:
function lap=laplaciano(f,v,h) % Calcula el Laplaciano de f en v, usando la formula simetrica central. % Se supone que f va de Rˆn a R, todos los vectores son columnas. % h es el parametro para calcular la derivada lap=0; n=length(v); for k=1:n e=zeros(n,1); e(k)=1; % k-esimo vector de la base canonica. lap=lap+(f(v+he)+f(v-he)-2*f(v))/hˆ2; end
Ejemplo 3.1.16 Escribe una f´ormula para calcular el rotacional de una funci´on f : R^3 → R^3 en el punto (0, 0 , 0) calculando ´unicamente la funci´on en 6 puntos. Respuesta: El rotacional de f viene dado por:
∇ × f =
∂fz ∂y −^
∂fy ∂z ∂fx ∂z −^
∂fz ∂x ∂fy ∂x −^
∂fx ∂y
Calculamos cada una de esas cantidades en (0, 0 , 0) con la f´ormula de la derivada central para derivadas parciales y para un h peque˜no de nuestra elecci´on:
∂fx ∂y
fx(0, h, 0) − fx(0, −h, 0) 2 h
∂fx ∂z
fx(0, 0 , h) − fx(0, 0 , −h) 2 h
∂fy ∂x
fy (h, 0 , 0) − fy (−h, 0 , 0) 2 h
∂fy ∂z
fy (0, 0 , h) − fy (0, 0 , −h) 2 h ∂fz ∂x
fz (h, 0 , 0) − fz (−h, 0 , 0) 2 h
∂fz ∂y
fz (0, h, 0) − fz (0, −h, 0) 2 h
Con ello, una f´ormula aproximada para el rotacional es:
∇ × f ≈
2 h
fz (0, h, 0) − fz (0, −h, 0) − (fy (0, 0 , h) − fy (0, 0 , −h))
fx(0, 0 , h) − fx(0, 0 , −h) − (fz (h, 0 , 0) − fz (−h, 0 , 0))
fy (h, 0 , 0) − fy (−h, 0 , 0) − (fx(0, h, 0) − fx(0, −h, 0))
N´otese que solo tenemos que evaluar la funci´on en 6 puntos para obtener toda esta informaci´on.
3.3. INTEGRALES DEFINIDAS 53
Exercise 3.
Write a program jacobianmatrix.m that on input f : Rn^ → Rm^ and x 0 ∈ Rn^ produces an approximation to the Jacobian matrix Jf (x 0 ) of f at x 0 using the five point formula for each direction. Use this program to compute the Jacobian matrix and the Jacobian determinant for the following functions.
f 1
x y
xy − e x 1+y^2 x
, f 2
x y z
x^2 + sin(yz) y − z^2
f 3
x y z
cos(x + y − z) sin(x^2 − y^2 + z^2 ) 1 1+x^2 +z^3
(^) , f 4
x y z
x y + x z + y + x z^2 − y^2 + xezy
Exercise 3.
Recall the chain rule for the Jacobian matrix of a composite function: if f : Rn^ → Rm^ and g : Rm^ → Rl^ then the Jacobian matrices satisfy: J(g ◦ f )(x 0 ) = Jg(f (x 0 ))Jf (x 0 ). (3.9)
Let f 1 , f 2 , f 3 , f 4 be as in 3.10. Compute the Jacobian matrix of f 1 ◦ f 1 using the left–hand term in (3.9) and your program of exercise 3.10. Now, compute it again but using the right–hand term of (3.10). Is the result exactly equal? The same question for f 1 ◦ f 2 , f 2 ◦ f 3 and f 4 ◦ f 3.
Exercise 3.
Using the data in exercise 2.14, approximate the value of the acceleration suffered by a person in free fall at different moments.
Exercise 3.
In exercise 5.8 you learnt how to produce realistic trajectories of cannon balls near to earth surface. Let us fix the parameters of the ball and air resistance in such a way that the mass of the ball is 10Kg and the modulus of the resistance force is 0, 01 ‖v‖ + 0, 001 ‖v‖^2 being v the speed. Produce a table which for a fixed angle of 45 degrees and different initial speeds from 100 to 500 meters per second describes:
The total distance (in the horizontal axis).
The final velocity of the ball in the moment of impact.
Then, compute the derivative of both quantities with respect to the speed of launch. Compute also the second derivative of the same function.
3.3. C´alculo aproximado de integrales definidas
Sea f : [a, b] → R una funci´on continua. De f supondremos ´unicamente que sabemos calcularla en valores concretos. Esto es, dado x 0 ∈ [a, b], podemos calcular f (x 0 ) (tal vez con un error dado) pero no tenemos a
nuestra disposici´on una f´ormula cerrada para f. Deseamos aproximar el valor de
∫ (^) b a f^ (x)^ dx. Alternativamente, puede sepamos la existencia de una funci´on pero ni siquiera podemos calcularla en los valores que deseemos, sino que disponemos ´unicamente de una tabla de valores. Nuestros m´etodos deben cubrir ambas necesidades.
54 CAP´ITULO 3. DERIVACI ON E INTEGRACI ´ ON NUM ´ ERICA APROXIMADA DE FUNCIONES´
3.3.1. M´etodos basados en la interpolaci´on
Partimos de que conocemos una tabla de datos y 1 = f (x 1 ),... , yk = f (xk) de la funci´on, sea porque esa es la ´unica informaci´on de la funci´on, sea porque podemos calcular la funci´on en cualquier punto que deseemos. Una familia de m´etodos para calcular integrales definidas se inspira en la siguiente idea.
Interpolar la funci´on en los valores dados, llamemos p(x) a la funci´on obtenida.
Calcular simb´olicamente la integral definida de p(x) entre a y b, siendo el valor calculado una aproximaci´on de la integral definida de f (x) entre a y b.
Afortunadamente, en muchos casos no es necesario hacer este proceso desde cero, sino que desarrollando las f´ormulas de antemano podemos obtener sencillas reglas de integraci´on y tener alguna informaci´on sobre los errores cometidos, sobre todo en el caso de que los x 1 ,... , xk sigan alg´un patr´on determinado. En las siguientes secciones desarrollamos algunas de estas reglas concretas en detalle.
Ejemplo 3.3.1 Usando los datos de la Figura 2.3 queremos aproximar el valor medio de la presi´on para tem- peraturas en el intervalo [2, 3 , 3 ,7]:
E(P (T )) =
2 , 3
P (T ) dT.
Para ello, aproximamos la integral de P (T ) (una funci´on desconocida) por la integral del polinomio p(x) que interpola a P (T ) en la tabla de datos que tenemos. Dicho polinomio fue calculado en el Ejemplo 2.1.3, y podemos calcular su integral de forma exacta obteniendo el valor:
∫ (^3) , 7
2 , 3
p(x) dx = 33, 3889 ,
con lo que concluimos:
E(P (T )) ≈
= 23,8492kP a.
Regla del punto medio (interpolaci´on polinomial)
Aplicamos el paradigma de esta secci´on con k = 1, x 1 = m el punto medio entre a y b. Consideramos por ello el polinomio interpolante en un punto:
p(x) = f (m) = f
a + b 2
Podemos entonces aproximar:
Teorema 3.3.2 (Regla del punto medio) Supongamos que f ∈ C^2 [a, b] y sea m = (a + b)/ 2 el punto medio entre a y b. Entonces, ∫ (^) b
a
f (x) dx = (b − a)f (m) +
(b − a)^3 24
f ′′(ζ),
para alg´un ζ ∈ (a, b).
56 CAP´ITULO 3. DERIVACI ON E INTEGRACI ´ ON NUM ´ ERICA APROXIMADA DE FUNCIONES´
El valor proporcionado por la regla del punto medio es
∫ (^) π/ 4
0
cos x dx ≈
π 4 cos
π 8
un error de 2 d´ecimas. Es un error considerable pero el trabajo que nos ha llevado es muy poco! Con m´etodos m´as avanzados conseguiremos mejores resultados.
Ejemplo 3.3.4 Aproximamos el valor de la integral de f (x) = e x x en el intervalo^ [1,^1 ,5]^ usando la regla del punto medio: (^) ∫ 1 , 5
1
ex x
dx ≈ (1, 5 − 1)
e^1 ,^25 1 , 25
El valor que da el software Matlab para la integral es de 1 , 4062 , con lo que nuestro error ronda las 3 d´ecimas.
Nota 3.3.5 Una demostraci´on independiente del Teorema 3.3.2 es como sigue: consideremos el polinomio in- terpolante p(x) de f (x) en m ± :
p(x) =
x − m − − 2 f (m − ) +
x − m + 2 f (m + ).
Es f´acil ver que lim→ 0 p(x) = f (m) + f ′(m)(x − m).
Se tiene por tanto usando el error de la interpolaci´on dado en el Teorema 2.1.6:
∫ (^) b
a
f (x) dx =
∫ (^) b
a
p(x) dx +
∫ (^) b
a
(x − m + )(x − m − ) 2
f ′′(ζ,x) dx.
La primera de las dos integrales en el t´ermino de la derecha tiende, en el l´ımite cuando → 0 , al valor (b − a)f (m). La segunda, en el mismo l´ımite, vale
∫ (^) b
a
(x − m)^2 2
f ′′(ζx) dx =
f ′′(ζ) 2
(b − a)^3 12
para alg´un ζ ∈ (a, b) usando el Teorema 1.2.2. El alumno avanzado puede notar aqu´ı el uso de un argumento de compacidad.
Regla del trapecio
Aplicamos el paradigma de esta secci´on con k = 2, x 1 = a, x 2 = b. El polinomio interpolante es en este caso:
p(x) = x − b a − b
f (a) + x − a b − a
f (b).
Teorema 3.3.6 (Regla del trapecio) Supongamos que f ∈ C^2 [a, b]. Entonces,
∫ (^) b
a
f (x) dx = (b − a)
f (a) + f (b) 2
(b − a)^3 12 f ′′(ζ),
para alg´un ζ ∈ (a, b).
3.3. INTEGRALES DEFINIDAS 57
Demostraci´on. Escribimos el polinomio interpolante con t´ermino de error:
f (x) = p(x) +
(x − a)(x − b) 2
f ′′(ζx),
donde ζx ∈ (a, b) depende de x. Integrando a ambos lados y utilizando el Teorema 1.2.2, ∫ (^) b
a
f (x) dx =(b − a)
f (a) + f (b) 2
∫ (^) b
a
(x − a)(x − b) 2
f ′′(ζx) dx
=(b − a)
f (a) + f (b) 2
- f ′′(ζ)
∫ (^) b
a
(x − a)(x − b) 2
dx,
para alg´un ζ ∈ (a, b). Calculando ex´actamente la ´ultima integral obtenemos el resultado deseado. Vemos una ilustraci´on gr´afica de este m´etodo en la Figura 3.4.
Figura 3.4: Ilustraci´on de la f´ormula del trapecio para aproximar integrales definidas. El ´area bajo la funci´on (esto es, la integral) se aproxima por el ´area del trapecio de la figura.
Ejemplo 3.3.7 Como en el Ejemplo 3.3.3, consideremos la integral de cos x entre 0 y π/ 4 , cuyo valor exacto es 0 , 7071. El valor proporcionado por la regla del trapecio es ∫ (^) π/ 4
0
cos x dx ≈
π 4
cos π 4 + cos 0 2
un error de 2 d´ecimas.
Ejemplo 3.3.8 Como en el Ejemplo 3.3.4, intentamos aproximar el valor de la integral de f (x) = e
x x en el intervalo [1, 1 ,5]. Ya hemos visto que el valor que da el software Matlab para la integral es de 1 , 4062. Usando la regla del trapecio: ∫ (^1) , 5
1
ex x
dx ≈ (1, 5 − 1)
e + e
1 , 5 1 , 5 2
con lo que nuestro error ronda las 2 cent´esimas.