Métodos numéricos PDF

Tema 4 Métodos numéricos Versión: 22 de septiembre de 2019 La mayor parte de las matemáticas estudiadas hasta ahora se han dedicado a desarrollar mét...

Tema 4 Métodos numéricos Versión: 22 de septiembre de 2019 La mayor parte de las matemáticas estudiadas hasta ahora se han dedicado a desarrollar métodos que nos proporcionen la solución exacta de un problema. Por ejemplo, calcular la solución de una ecuación del tipo f (x) = 0 realizando operaciones elementales sobre la misma para conseguir despejar la incógnita x. Z b O bien, calcular el valor de una integral definida f (x) dx calculando una primitiva de f (x) y luego aplicando a la Fórmula de Barrow. Desgraciadamente, en la gran mayoría de los casos que se presentan en la práctica, estos métodos no son de aplicación. Ello puede deberse a que el método para calcular la solución exacta sea muy complicado, a que no se conozca un método adecuado, o incluso a que no exista un método que nos permita, mediante cálculos elementales, encontrar la solución. En estos casos es necesario recurrir a métodos numéricos, denominados así porque, usualmente, consisten en realizar una sucesión más o menos larga de operaciones numéricas (normalmente mediante la ayuda de un ordenador), al cabo de las cuales encontramos un valor numérico que, si bien no es la solución exacta del problema, se le parece mucho, es decir, aproxima la solución buscada con una precisión razonablemente buena. 4.1 Resolución numérica de ecuaciones Uno de los problemas que más se presenta en matemáticas es el de calcular la solución de una ecuación. En algunas (pocas) ocasiones, esto puede hacerse por métodos analíticos, es decir, se puede “despejar” la incógnita para encontrar el o los valores que resuelven la ecuación. En la gran mayoría de las ocasiones con algún interés práctico esto no es posible y es necesario recurrir a un método numérico que, con la ayuda de un ordenador, nos permita calcular un valor aproximado de la solución. 4.1.1 Teoremas del Valor Intermedio y de Bolzano Cuando se plantea el problema de calcular la solución de una ecuación como f (x) = 0 existen dos cuestiones previas que conviene analizar: ¿Tiene solución esta ecuación? ¿Dónde está (aunque sea más o menos) la solución? En los casos en que la solución se puede calcular exactamente por métodos elementales, se tiene la respuesta a ambas preguntas: existe, puesto que la hemos encontrado y sabemos dónde está, puesto que sabemos exactamente su valor. 159 4. Métodos numéricos 160 En muchos de los otros casos, la respuesta a estas preguntas se obtiene con ayuda de los siguientes teoremas. Teorema del Valor Intermedio Una función continua en un intervalo [a, b] toma todos los valores comprendidos entre f (a) y f (b). Figura 4.1: Teorema del Valor Intermedio: como se Figura 4.2: Teorema del Valor Intermedio: además puede observar, la función toma todos los valores de todos los valores comprendidos entre f (a) y f (b) comprendidos entre f (a) y f (b). la función f puede tomar otros valores. Teorema de Bolzano Sea f una función continua en un intervalo [a, b] y tal que f (a) y f (b) tienen signos opuestos (es decir f (a)f (b) < 0). Entonces existe c 2 (a, b) tal que f (c) = 0. Figura 4.3: Teorema de Bolzano: como se puede observar, la fun- ción tiene signos opuestos en a y b (f (a) < 0 y f (b) > 0). En concecuencia, toma el valor 0 en algún punto del intervalo (a, b) (de hecho lo toma en tres puntos). Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 161 Ejemplo 4.1 Utilizando el Teorema de Bolzano, probar que la ecuación x = 2 x tiene al menos una solución real. En primer lugar, hay que escribir la ecuación en la forma f (x) = 0 y, luego, encontrar un intervalo [a, b] en el cual se verifiquen las hipótesis del Teorema, para así poder concluir que existe algún punto en el intervalo en el que la función se anula, es decir, alguna solución de la ecuación. Se tiene: x x x=2 () f (x) = x 2 =0 Esta función está definida y es continua en todo R. Es fácil ver que f (0) = 0 20 = 1 < 0. Por otro lado, teniendo en cuenta que cuando x tiende a +1, lı́m 2 x = 0, tampoco es difícil comprender que, para x x!+1 suficientemente grande, x será mayor que 2 x y por tanto x 2 x será positivo. Por ejemplo: 1 1 1 f (1) = 1 2 = >0=1 2 2 En consecuencia, f verifica las hipótesis del Teorema de Bolzano en el intervalo [0, 1]: es continua y f (0) y f (1) tienen signos opuestos. Luego podemos afirmar que f (x) tiene al menos un cero en el intervalo (0, 1). O, lo que es lo mismo, que la ecuación x = 2 x tiene al menos una solución en dicho intervalo. Ejemplo 4.2 Utilizando el Teorema de Bolzano, probar que la ecuación x4 = 1 + 3e x tiene al menos una raíz real. Razonando como en el ejercicio anterior, se tiene x4 = 1 + 3e x () f (x) = x4 1 3e x =0 La función f (x) está definida y es continua en todo R. Se tiene, por ejemplo, f (0) = 1 3 = 4 < 0. Por otro lado, igual que en el ejemplo anterior, x4 1 tiende a +1 cuando x ! +1 mientras que lı́m 3e x = 0, y no resulta difícil comprender que, para x suficientemente grande, x4 1 será mayor que x!+1 3e y por tanto x4 1 3e x será positivo. x Por ejemplo, recordando que e x < 1 8x > 0 y, en consecuencia, que 3e x < 3 8x > 0, se tiene f (2) = 24 1 3e 2 = 15 3e 2 > 12 > 0 Luego, por el Teorema de Bolzano, f (x) tiene, al menos, un cero en el intervalo (0, 2), es decir, la ecuación dada tiene, al menos, una raíz en dicho intervalo. 4.1.2 Resolución numérica de ecuaciones: método de bisección Se presenta aquí un método sencillo, basado directamente en el Teorema de Bolzano, que permite, en determi- nadas circunstancias, calcular la solución de una ecuación. Hay que comenzar por decir que cualquier ecuación en una variable se puede siempre escribir (y no de manera única) en la forma de una equivalente (es decir, que tiene las mismas soluciones) pero con segundo miembro nulo f (x) = 0 Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 162 Ejemplo 4.3 La ecuación x = 2 se puede también escribir x 2 x = 0. x 1 También se tiene x = 2 x , x = x , x 2x = 1, luego también se puede escribir x 2x 1 = 0. 2 Dada f : [a, b] ⇢ R 7! R, continua, se plantea el problema de encontrar una solución (también llamada raíz) de la ecuación f (x) = 0. Desde el punto de vista geométrico, esto significa encontrar, en [a, b], un punto de corte de la gráfica de la función y = f (x) con el eje de abscisas (ver la Figura 4.4). a ↵ b Figura 4.4: La gráfica de y = f (x) corta al eje de abscisas en un punto ↵ 2 [a, b], lo que significa que ↵ es una solución de la ecuación f (x) = 0. Los métodos de aproximación de raices de ecuaciones necesitan conocer, o bien un intervalo que contenga sólo una raíz, o bien un punto inicial que esté suficientemente cerca de ella. Por tanto, como paso previo a la aplicación de un método de aproximación, es necesario localizar la raíz, es decir encontrar un intervalo que la contenga y separar la raíz, es decir encontrar un intervalo que sólo contenga dicha raíz. Esto se hace por métodos analíticos, gráficos y, en algunos casos, empíricos. Ejemplo 4.4 0.5 0 −0.5 ↵ −1 −1.5 −2 −2.5 −3 −1 −0.5 0 0.5 1 La función y = x 2 x = f (x) está representada en la Figura para x 2 [ 1, 1]. Se observa que hay un único punto ↵ 2 [0, 1] en que la curva corta al eje OX, es decir, que hay una única raiz de x 2 x = 0 en [0, 1]. ( f (0) = 0 20 = 1 < 0, 1 1 f (1) = 1 2 1 = 1 = > 0. 2 2 Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 163 Los métodos para aproximar raíces de ecuaciones son, en general iterativos, es decir consisten en construir una sucesión de valores x1 , x2 , x3 , x4... mediante una relación de recurrencia, esto es, se calcula cada uno de ellos a partir del anterior: x1 ! x2 ! x3 ! x4 , etc. Cuando la sucesión de valores x1 , x2 , x3... tiende hacia la raíz ↵ de f (es decir, se acerca cada vez más a ella, tanto como se quiera: lı́m xn = ↵), se dice que el método iterativo es convergente. n!1 Método de bisección Sin mucha precisión, el método de bisección consiste en lo siguiente: 1. Subdividir en dos partes el intervalo en que se sabe que la función cambia de signo y tiene una sola raíz. 2. Averiguar, utilizando el Teorema de Bolzano, en cual de las dos mitades se encuentra la raiz y descartar la otra mitad del intervalo. 3. Reiniciar este proceso con el subintervalo elegido. 4. Continuar con este proceso hasta que el subintervalo elegido tenga una longitud lo suficientemente pequeña como para que cualquiera de sus puntos sea una aproximación aceptable de la solución. La elección óptima como aproximación es, entonces, el punto medio del subintervalo. a ↵ b a ↵ b a ↵ b x1 x2 x3 Figura 4.5: Tres etapas del método de dicotomía. En cada iteración se descarta la mitad del intervalo que no contiene a la raíz (en la que f no cambia de signo). El intervalo donde se encuentra la raíz es cada vez más pequeño y, su punto medio se acerca cada vez más a la solución buscada. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 164 Ejemplo 4.5 Utilizando el método de dicotomía, aproximar la solución de la ecuación x 2 x = 0 en el intervalo [0, 1] Sea f (x) = x 2 x. Intervalo Punto medio 0+1 [0, 1] f (0) < 0 f (1) > 0 x0 = = 0.5 2 0.5 + 1 [0.5, 1] f (0.5) < 0 x1 = = 0.75 2 0.5 + 0.75 [0.5, 0.75] f (0.75) > 0 x2 = = 0.625 2 0.625 + 0.75 [0.625, 0.75] f (0.625) < 0 x3 = = 0.6875 2 0.625 + 0.6875 [0.625, 0.6875] f (0.6875) > 0 x4 = = 0.65625 2... Por lo que una aproximación de la solución es ↵ ⇡ 0.65625, obtenida aplicando el proceso de subdivisión 4 veces y eligiendo como aproximación el punto medio del último subintervalo. Obsérvese que, si se elige como aproximación x0 , el error máximo que se comete es la mitad de la longitud b a del intervalo inicial e0 =. Si se elige como aproximación x1 , el error máximo es la mitad del anterior 2 e0 b a b a e1 = =. Reiterando este razonamiento, si se elige como aproximación xn , el error máximo es en = n+1. 2 22 2 Esto permite saber, a priori, cuantas iteraciones hay que realizar para conseguir una aproximación con un error tan pequeño como se quiera. En efecto, si en el intervalo [a, b] hay una solución ↵, ¿qué número n de veces hay que aplicar el proceso de subdivisión para conseguir que el error cometido no sea mayor que una cantidad dada "? Se ha visto que, si se aplica n veces, el error máximo que se comete tomando xn como aproximación es b a en = 2n+1 En consecuencia habrá que elegir n de forma que se tenga ✓ ◆ b a ✓ ◆ ln b a b a b a " n+1 = = ⇡ 6.64 () n > 6.64 1 = 5.64 ln(2) ln(2) ln(2) Luego hay que realizar 6 iteraciones. Si se hicieran estas 6 iteraciones se obtendría como aproximación x6 = 0.6484, cuyas dos primeras cifras decimales son exactas, y por lo tanto, el error de aproximación es menor que 0.01. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 165 Ejemplo 4.7 Utilizando el método de dicotomía, aproximar la solución de la ecuación del Ejercicio 4.2, x4 = 1 + 3e x , en el intervalo [0, 2] con un error menor que 0.05 16 14 12 10 8 6 4 2 ↵ 0 −2 −4 0 0.5 1 1.5 2 Sea f (x) = x4 1 3e x. Como se puede observar en la figura, f tiene una única raíz en [0, 2]. Puesto que se desea un error menor que 0.05, habrá que tomar ✓ ◆ ✓ ◆ b a 2 ln ln " 0.05 ln(40) n+1> = = ⇡ 5.32 () n > 4.32 ln(2) ln(2) ln(2) Luego hay elegir n = 5 (es decir, elegir como aproximación x5 ). intervalo pto. medio error 0+2 [0, 2] f (0) < 0 f (2) > 0 x0 = =1 1 2 1+2 [1, 2] f (1) < 0 x1 = = 1.5 0.5 2 1 + 1.5 [1, 1.5] f (1.5) > 0 x2 = = 1.25 0.25 2 1 + 1.25 [1, 1.25] f (1.25) > 0 x3 = = 1.125 0.125 2 1.125 + 1.25 [1.125, 1.25] f (1.125) < 0 x4 = = 1.1875 0.0625 2 1.125 + 1.875 [1.125, 1.1875] f (1.1875) > 0 x5 = = 1.15625 0.03125 2 Por lo que una aproximación de la solución es ↵ ⇡ 1.15625 con un error menor o igual que 0.05 4.1.3 Resolución numérica de ecuaciones: método de Newton El método de bisección, presentado en la sección anterior, sólo hace uso de los valores que toma la función f cuyas raíces se quieren calcular. En esta sección se presenta un método que utiliza además los valores que toma la derivada de f. Naturalmente, esto requiere que f sea derivable. Sea, pues, f : [a, b] ⇢ R ! R una función continua, derivable y con derivada continua. Se supone que la ecuación f (x) = 0 tiene en el intervalo (a, b) una única solución ↵, que no se conoce y se desea aproximar: f (↵) = 0, ↵ 2 (a, b) Se recuerda que ↵ es un punto de corte de la gráfica de y = f (x) con el eje OX. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 166 a b ↵ La idea del método de Newton consiste en sustituir, en determinados puntos, la gráfica de la función por la de su recta tangente en dichos puntos. Se comienza eligiendo un punto inicial x0 2 [a, b], que debe estar cerca de la solución ↵ que se quiere aproximar. La ecuación de la recta tangente a y = f (x) en el punto (x0 , f (x0 )) es (ver Figura 4.6) y = f (x0 ) + f 0 (x0 )(x x0 ) (x 0 , f (x 0 )) (x 0 , f (x 0 )) a b a b ↵ x0 ↵ x1 x0 Figura 4.6: La recta tangente a la curva y = f (x) en Figura 4.7: La recta tangente a la curva y = f (x) en el punto (x0 , f (x0 )) tiene de ecuación el punto (x0 , f (x0 )) corta al eje OX en y = f (x0 ) + f 0 (x0 )(x x0 ). f (x0 ) x1 = x0. f 0 (x0 ) Esta recta corta al eje OX en el punto de abscisa f (x0 ) x1 = x0 f 0 (x0 ) Parece claro que el punto x1 está más cerca de ↵ que el punto inicial x0. Se repite ahora el proceso anterior, pero comenzando en el punto x1. El método de Newton consiste en reiterar este proceso, partiendo cada vez del punto calculado en la etapa anterior. Esto proporcionará puntos cada vez más cercanos a la solución ↵. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 167 (x 0 , f (x 0 )) (x 0 , f (x 0 )) (x 1 , f (x 1 )) (x 1 , f (x 1 )) a b a b ↵ x1 x0 ↵ x2 x1 x0 Figura 4.8: La recta tangente a la curva y = f (x) en Figura 4.9: La recta tangente a la curva y = f (x) en el punto (x1 , f (x1 )) tiene de ecuación el punto (x1 , f (x1 )) corta al eje OX en y = f (x1 ) + f 0 (x1 )(x x1 ). f (x1 ) x2 = x1. f 0 (x1 ) Método de Newton Consiste en lo siguiente: 1. Elegir un punto x0 que esté cerca de la solución. 2. Calcular sucesivamente los puntos f (xn ) xn+1 = xn , para n = 0, 1, 2,... f 0 (xn ) hasta conseguir una aproximación lo suficientemente buena de ↵. Observaciones: 1. En la descripción anterior hay dos indefiniciones claras: a) ¿Cómo se elige un punto x0 que esté cerca de la solución? No hay una respuesta general a esta pregunta. Puede que se conozca, por ejemplo, por razones empíricas o por análisis previo. Si no, una posibilidad es utilizar previamente el método de bisección y comenzar el método de Newton en la solución proporcionada por aquél. b) ¿Cómo se sabe si una aproximación es lo suficientemente buena? En la práctica, lo que se suele hacer cuando se utiliza este método con un ordenador, es detenerse cuando dos aproximaciones consecutivas están muy cercanas: |xn+1 xn | < una cantidad muy pequeña previamente fijada, por ejemplo 10 4 2. Como se ha visto, en el método de Newton hay que dividir por el valor de la derivada de f en determinados puntos, que están cercanos a la solución. Naturalmente, es imprescindible, pues, que la derivada f 0 no se anule cerca de la solución. 3. Este método utiliza mucha más información sobre la función f que el método de bisección, que se vió en el Tema 3, ya que hace uso de la derivada. Es por ello lógico que sea mejor, es decir más rápido en llegar a la solución. De hecho es mucho más rápido. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 168 Ejemplo 4.8 Determinar el número de soluciones en R de la ecuación siguiente y utilizar el método de Newton para aproximar la mayor de ellas. ex +x 2 = 0 a) Denotemos f (x) = ex +x 2. Sabemos que f es derivable en R y f 0 (x) = ex +1 > 0 8x 2 R lo cual significa que f es creciente en R. 10 También se tiene 8 lı́m f (x) = +1 y lı́m f (x) = 1 6 x!+1 x! 1 y = ex + x 2 4 Gráficamente deducimos que f sólo tiene una raíz, es decir, 2 la ecuación f (x) = 0 tiene una única solución ↵ 2 R. 0 Como f (0) = 1 < 0 y f (1) = e 1 > 0, por el Teorema de Bolzano se tiene que la raíz está en el intervalo (0, 1). −2 −4 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 b) Utilizamos ahora el método de Newton para aproximar ↵. Tomamos como primer punto x0 = 0. Se tiene: f (x0 ) e0 + 0 2 1 x1 = x0 =0 = = 0.5 f 0 (x0 ) e0 + 1 2 Partiendo de x1 , calculamos f (x1 ) e0.5 + 0.5 2 x2 = x1 = 0.5 ⇡ 0.44385167 f 0 (x1 ) e0.5 + 1 Repetimos el proceso y calculamos f (x2 ) x3 = x2 ⇡ 0.44285470 f 0 (x2 ) Repetimos el proceso una vez más y obtenemos f (x3 ) x4 = x3 ⇡ 0.44285440 f 0 (x3 ) Observamos que las 6 primeras cifras decimales de las dos últimas aproximaciones son iguales: 0.442854, de manera que se tiene: |x4 x3 | = 0.00000030 = 3 ⇥ 10 7 < 10 6 Tomamos, pues x4 = 0.44285440 como aproximación de la solución. Observación: Hacer estos cálculos a mano no es sencillo. Pero sí lo es hacerlos con una hoja de cálculo EXCEL. Es interesante hacerlo como ejercicio. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 169 Ejemplo 4.9 Utilizando el método de Newton, aproximar la solución de la ecuación x 2 x = 0 en el intervalo [0, 1]. Denotemos f (x) = x 2 x. Sabemos que f es derivable en R y f 0 (x) = 1 + ln(2) 2 x Utilizamos ahora el método de Newton para aproximar la solución de la ecuación. Tomamos como primer punto x0 = 0. Se tiene: f (x0 ) 0 20 x1 = x0 0 =0 ⇡ 0.590616109 f (x0 ) 1 + ln(2) 20 Partiendo de x1 , calculamos f (x1 ) x2 = x1 ⇡ 0.640909617 f 0 (x1 ) Repetimos el proceso y calculamos f (x2 ) x3 = x2 ⇡ 0.641185736 f 0 (x2 ) Repetimos el proceso una vez más y obtenemos f (x3 ) x4 = x3 ⇡ 0.641185744 f 0 (x3 ) Observamos que las 7 primeras cifras decimales de las dos últimas aproximaciones son iguales: 0.6411857. De hecho esto indica, en general, que dichas 7 primeras cifras son exactas (en este caso, en concreto, todas las cifras de x4 son exactas). Se tiene: 9 8 |x4 x3 | = 0.000000008 = 8 ⇥ 10 < 10 Tomamos, pues x4 = 0.6411857 como aproximación de la solución. Observación: Esta misma ecuación fue resuelta, en el Ejercicio 4.5, por el método de bisección, encontrándose allí la aproximación 0.65625 tras 4 iteraciones. Esta aproximación sólo tiene una cifra decimal exacta: 0.6. Con el método de Newton hemos encontrado una aproximación con 9 cifras decimales exactas en 4 iteraciones. Resulta obvio, pues, que este método es (mucho) más rápido que el de bisección (de hecho es el más rápido). Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 170 4.2 Nociones de integración numérica Como se ha visto antes, si se conoce una primitiva F de la función f , se puede calcular el valor de la integral definida mediante la Regla de Barrow: Z b f (x) dx = F (b) F (a). a En la mayoría de los casos, sin embargo, no se puede utilizar esta fórmula, ya que no se conoce dicha primitiva. Es posible, por ejemplo, que no se conozca la expresión matemática de la función f , sino sólo sus valores en determinados puntos, recogidos de un experimento. Pero también hay funciones (de apariencia sencilla) para las que se puede demostrar que no tienen ninguna primitiva que pueda escribirse en términos de funciones 2 elementales (por ejemplo e x ) La integración numérica es una herramienta de las matemáticas que proporciona fórmulas y técnicas para calcular aproximaciones de integrales definidas. Gracias a ella se pueden calcular, bien es cierto que de forma aproximada, valores de integrales definidas que no pueden calcularse analíticamente y, sobre todo, se puede realizar ese cálculo en un ordenador. Z b La idea básica para aproximar el valor de f (x) dx sin utilizar una primitiva de f ya se expuso en la sección 3.6: a calcular la suma de las áreas de los rectángulos que “recubren” el área. y y y y=f(x) y=f(x) y=f(x) x a b a b x a b x Z b Figura 4.10: La integral definida f (x) dx , que es el valor del área bajo la curva sombreada en la primera a figura, se puede aproximar por el resultado de sumar las áreas de los rectángulos. Como resulta evidente, se comete un error, ya que se desprecian –en este caso– las áreas de las pequeñas zonas triangulares comprendidas entre la curva y los rectángulos. En el caso particular de la función representada en las figuras, el valor de la aproximación es menor que el valor exacto. Pero en otros casos puede ser mayor; véase, por ejemplo, la figura siguiente. y a b x Figura 4.11: En este caso,la suma de las áreas de los rectángulos proporciona un valor mayor que el valor exacto, pero igualmente es una aproximación. Como también resulta evidente, y se puede demostrar matemáticamente, el error que se comete es más pequeño (en valor absoluto, es decir, sin tener en cuenta el signo del mismo) cuanto más “estrechos” sean los rectángulos, es decir, cuanto mayor cantidad de ellos se usen. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 171 ¿Cómo se calcula la suma de las áreas de los rectángulos? Se supone que se usan 5 rectángulos, como en la Figura 4.12 y se denotan x1 = a, x2 , x3 , x4 , x5 y x6 = b los puntos que determinan los 5 subintervalos. Se supone también, para hacer las cosas más fáciles, que estos puntos están regularmente espaciados, es decir, que la distancia entre cada dos puntos consecutivos, que se denota h, es siempre la misma. El área de los distintos rectángulos es (recordando área = base⇥altura): Area(R1 ) = Longitud del segmento [x1 , x2 ] ⇥ Altura del rectángulo = (x2 x1 ) ⇥ f (x1 ) = h f (x1 ) Area(R2 ) = Longitud del segmento [x2 , x3 ] ⇥ Altura del rectángulo = (x3 x2 ) ⇥ f (x2 ) = h f (x2 ) etc. Sumando todas se tiene: Area(R1 ) + · · · + Area(R5 ) = hf (x1 ) + hf (x2 ) + hf (x3 ) + hf (x4 ) + hf (x5 ) ⇣ ⌘ = h f (x1 ) + f (x2 ) + f (x3 ) + f (x4 ) + f (x5 ) y esta última expresión proporciona una aproximación (es verdad que no muy buena, de momento) del valor de la integral: Z b ⇣ ⌘ f (x) dx ⇡ h f (x1 ) + f (x2 ) + f (x3 ) + f (x4 ) + f (x5 ) a Observamos ahora que, puesto que hay 5 subintervalos de igual longitud, debe ser Longitud del intervalo [a, b] b a h= = 5 5 luego, la fórmula anterior quedaría Z b b a⇣ ⌘ f (x) dx ⇡ f (x1 ) + f (x2 ) + f (x3 ) + f (x4 ) + f (x5 ) a 5 y y=f(x) f(x5) f(x4) f(x3) f(x2) f(x1) R R 5 R 4 R 3 R 2 1 h a= x x2 x3 x4 x5 b=x6 x 1 Figura 4.12: La altura del rectángulo de base [x1 , x2 ] es f (x1 ), el valor de f en x1 ; la del rectángulo de base [x2 , x3 ] es f (x2 ); etc. Si, en lugar de 5, tuviéramos 6 subintervalos, entonces tendríamos 7 puntos: x1 = a, x2 , x3 , x4 , x5 , x6 y x7 = b y la aproximación se escribiría: Z b b a⇣ ⌘ f (x) dx ⇡ f (x1 ) + f (x2 ) + f (x3 ) + f (x4 ) + f (x5 ) + f (x6 ) a 6 (obsérvese que el último punto x7 no se utiliza en esta expresión). Si el número de subintervalos utilizados fuera muy grande, por ejemplo, 100 (es decir, 101 puntos), se podría escribir Z b b a⇣ ⌘ f (x) dx ⇡ f (x1 ) + f (x2 ) + · · · + f (x100 ) a 100 Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 172 Es preferible y más usual, sin embargo, utilizar la expresión siguiente Z 100 b b a X f (x) dx ⇡ f (xi ) a 100 i=1 P El símbolo (letra griega sigma mayúscula) es muy utilizado en matemáticas: se denomina “sumatorio” y sirve para escribir de forma escueta una suma con un número muy grande o indeterminado de sumandos. 100 X La expresión f (xi ) se lee : suma de f (xi ) desde i = 1 hasta i = 100. i=1 Ya podemos, pues, escribir de forma general la aproximación de la integral para un número indeterminado de subintervalos. Fórmula de los rectángulos Sea f una función continua en [a, b] y sean x1 = a, x2 , x3 ,... , xn+1 = b, n + 1 puntos que definen una partición b a del intervalo [a, b] en n subintervalos, todos de la misma longitud h =. n Entonces la integral definida de f entre a y b se puede aproximar por Z n b b a X f (x) dx ⇡ f (xi ) a n i=1 En la deducción de esta fórmula se ha aproximado el área bajo la curva en cada subintervalo por el área del rectángulo con la misma base y altura igual al valor de la función en el extremo inferior del subintervalo, como en la Figura 4.13. Pero también se podría haber utilizado el valor de la función en el extremo superior, como se ve en la Figura 4.14. y y x1 x2 x x1 x2 x Figura 4.13: Se toma como altura del rectán- Figura 4.14: Se toma como altura del rectán- gulo el valor de f en el extremo inferior, x1. gulo el valor de f en el extremo superior, x2. Así se obtendría una variante de la Fórmula de los Rectángulos. Ambas fórmulas dan resultados similares desde el punto de vista del error que se comete en la aproximación. Fórmula de los rectángulos (variante) Sea f una función continua en [a, b] y sean x1 = a, x2 , x3 ,... , xn+1 = b, n + 1 puntos que definen una partición b a del intervalo [a, b] en n subintervalos, todos de la misma longitud h =. n Entonces la integral definida de f entre a y b se puede aproximar por Z n n+1 b b a X b a X f (x) dx ⇡ f (xi+1 ) = f (xi ) a n i=1 n i=2 Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 173 Otra posibilidad, es tomar como altura del rectángulo el valor de la función en el punto medio del subintervalo, como se muestra en la Figura 4.15 y x1 x 1+ x 2 x2 x 2 Figura 4.15: En la Fórmula del punto medio, se aproxima el área bajo la curva por el área del rectángulo de altura igual al valor de la función en el punto medio del subintervalo. Fórmula del punto medio Sea f una función continua en [a, b] y sean x1 = a, x2 , x3 ,... , xn+1 = b, n + 1 puntos que definen una partición b a del intervalo [a, b] en n subintervalos, todos de la misma longitud h =. n Entonces la integral definida de f entre a y b se puede aproximar por Z n ✓ ◆ b b a X xi + xi+1 f (x) dx ⇡ f a n i=1 2 Esta fórmula es de orden 1. y y y=f(x) y=f(x) a= x x2 x3 x4 x5 b=x6 x a= x x2 x3 x4 x5 b=x6 x 1 1 Figura 4.16: Fórmula de los rectángulos to- Figura 4.17: En la Fórmula del punto medio mando como altura el valor de f en el extre- elige como altura de los rectángulos en valor mos superior de cada subintervalo. de la función los puntos medios de cada subin- tervalo. Orden de una fórmula de integración numérica Se dice que una fórmula de integración es de orden k cuando es exacta para polinomios de grado k, es decir, que cuando el integrando es un polinomio de grado k, la fórmula proporciona el valor exacto de la integral. El orden de una fórmula de integración numérica nos da una medida de su bondad. La Fórmula de los rectángulos es de orden 0. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 174 Ejemplo 4.10 Z 1 x2 Aproximar el valor de la integral definida e dx utilizando la fórmula de los rectángulos 1 con 8 subintervalos. Se construye una partición de [ 1, 1] en 8 subintervalos, de forma que y 1 ( 1) 2 1 h= = = = 0.25 8 8 4 y=e−x 2 y los puntos del soporte de la partición son: x1 = 1 = 1 x6 = 1 + 5h=0.25 x2 = 1+h = 0.75 x7 = 1 + 6h=0.5 x3 = 1 + 2h= 0.5 x8 = 1 + 7h=0.75 a=x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 b=x 9 x x4 = 1 + 3h= 0.25 x9 = 1 + 8h=1 x5 = 1 + 4h= 0 Según la Fórmula de los Rectángulos anterior: Z 1 8 X x2 x2i e dx ⇡ h e 1 i=1 Con ayuda de una calculadora, se tiene: Z 1 ⇣ ⌘ 2 e x dx ⇡ 0.25 0.3679 + 0.5698 + 0.7788 + 0.9394 + 1 + 0.9394 + 0.7788 + 0.5698 = 1.4860 1 Hay que insistir en que el valor calculado es sólo una aproximación del valor de la integral definida. Otra posibilidad es aproximar el área bajo la curva en cada subintervalo por el área del trapecio que se muestra en la Figura 4.18. y y y=f(x) f (x 2 ) f (x 1 ) h x1 x2 x a= x x2 x3 x4 x5 b=x6 x 1 Figura 4.18: En el subintervalo [x1 , x2 ], por Figura 4.19: En la Fórmula de los trapecios, ejemplo, el área bajo la curva se aproxima por se aproxima el valor de la integral definida por el área del trapecio, que tiene una base de lon- la suma de las áreas de los trapecios. gitud f (x1 ), otra base de longitud f (x2 ), y al- tura h = x2 x1. suma de las bases Recordando que el área de un trapecio es = ⇥ altura, se tiene que el área del trapecio 2 de la Figura 4.18 es f (x1 ) + f (x2 ) h 2 y que la suma de las áreas de todos los de la Figura 4.19, es decir la aproximación de la integral, es Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 175 Z b f (x1 ) + f (x2 ) f (x2 ) + f (x3 ) f (x5 ) + f (x6 ) f (x) dx ⇡ h+ h + ··· + h a 2 2 2 h⇣ ⌘ = f (x1 ) + f (x2 ) + f (x2 ) + f (x3 ) + · · · + f (x5 ) + f (x6 ) 2 b a⇣ ⌘ = f (x1 ) + 2f (x2 ) + 2f (x3 ) + 2f (x4 ) + 2f (x5 ) + f (x6 ) 2⇥5 Obsérvese que, en esta suma, el valor de f en los extremos (x1 = a y x6 = b) aparece una sola vez, mientras que el valor en los puntos internos (x2 , x3 , x4 y x5 ) aparece dos veces. Generalizando esto al caso general, con un número indeterminado de subintervalos, se tiene: Fórmula de los trapecios Sea f una función continua en [a, b] y sean x1 = a, x2 , x3 ,... , xn+1 = b, n + 1 puntos que definen una partición b a del intervalo [a, b] en n subintervalos, todos de la misma longitud h =. n Entonces la integral definida de f entre a y b se puede aproximar por Z b n ! b a X f (x) dx ⇡ f (a) + 2 f (xi ) + f (b) a 2n i=2 Esta fórmula es de orden 1. Ejemplo 4.11 Z 1 2 Aproximar el valor de la integral definida sen(ex ) dx utilizando la fórmula de los trapecios 0 con 5 subintervalos. Se considera una partición de [0, 1] en 5 subintervalos, de forma que y 1 h = = 0.2 5 y los puntos del soporte de la partición son: x1 = 0 x21 = 0 x2 = 0.2 x22 = 0.04 x3 = 0.4 x23 = 0.16 a=x1 x2 x3 x4 x5 b=x6 x x4 = 0.6 x24 = 0.36 x5 = 0.8 x25 = 0.64 x6 = 1 x26 = 1 La Fórmula de los trapecios anterior: Z " 5 # 1 h X x2 0 x2i 1 sen(e ) dx ⇡ sen(e ) + 2 sen(e ) + sen(e ) 0 2 i=2 ⇥ ⇤ = 0.1 sen(e0 ) + 2 sen(e0.04 ) + 2 sen(e0.16 ) + sen(e0.36 ) + 2 sen(e0.64 ) + sen(e1 ) Se tiene: Z 1 h i 2 sen ex dx ⇡ 0.1 0.8415 + 2 0.8628 + 0.9221 + 0.9906 + 0.9474 + 0.4108 1 = 0.8698 Hay que insistir en que el valor calculado es sólo una aproximación del valor de la integral definida. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 176 4.3 Interpolación y ajuste de datos En ciencias experimentales con frecuencia es necesario trabajar con conjuntos discretos de valores de alguna magnitud que depende de otra variable. Pueden proceder de muestreos, de experimentos o incluso de cálculos numéricos previos. Por ejemplo, se puede disponer de unos valores obtenidos experimentalmente sobre el número de individuos de una determinada especie de peces en un lago, obtenidos en distintos momentos a lo largo de un año. En ocasiones, para utilizar estos valores en cálculos posteriores es preciso «darles forma» de función, es decir: es preciso disponer de una función dada por una expresión matemática que «coincida» con dichos valores. Por ejemplo, se puede querer saber el número de peces que había en el lago en un momento intermedio para el que no se dispone de datos. Existen básicamente dos enfoques para conseguir esto: Interpolación es el proceso de determinar una función que tome exactamente los valores dados para los valores adecuados de la variable independiente, es decir que pase exactamente por unos puntos dados. Por ejemplo, determinar un polinomio de grado 4 que pase por 5 puntos dados, como en la figura de la derecha. Ajuste de datos es el proceso de determinar la función, de un tipo determinado, que mejor se aproxime a los datos («mejor se ajuste»), es decir tal que la distancia a los puntos (medida de alguna manera) sea lo menor posible. Esta función no pasará necesariamente por los puntos dados. Por ejemplo, determinar un polinomio de grado 1 que aproxime lo mejor posible unos datos, como se muestra en la figura adjunta. Cuando se trata de interpolar por un polinomio de un determinado grado, se habla de interpolación polinómica. Interpolación lineal Es sabido que por dos puntos dados del plano, (x1 , y1 ) y (x2 , y2 ), con x1 6= x2 , pasa una sola línea recta. Sea y = ax + b su ecuación. Se trata de determinar los valores que deben tener a y b para que, efectivamente, esa recta pase por esos puntos. Para ello se tiene que verificar: ⇢ y1 = ax1 + b y2 = ax2 + b La solución de este sistema lineal de dos ecuaciones con dos incógnitas proporciona los valores adecuados de los coeficientes a y b. y (x 2 , y 2 ) x (x 1 , y 1 ) Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 177 Interpolación cuadrática En general, por tres puntos dados del plano, (x1 , y1 ), (x2 , y2 ) y (x3 , y3 ), pasa una única parábola (polinomio de grado 2). Sea y = ax2 + bx + c su ecuación. Para calcular los valores adecuados de los coeficientes hay que resolver el sistema lineal de ecuaciones 8 2 2 32 3 2 3 < y1 = ax21 + bx1 + c x1 x1 1 a y1 y2 = ax22 + bx2 + c que, en forma matricial es 4 x22 x2 1 5 4 b 5 = 4 y2 5 : 2 y3 = ax3 + bx3 + c x23 x3 1 c y3 y su solución (única) proporciona los coeficientes que determinan la función interpolante. y (x 3 , y 3 ) (x 1 , y 1 ) x (x 2 , y 2 ) Interpolación polinómica global En general, dados N puntos (xk , yk ), k = 1,... , N , con xk todos distintos, existe un único polinomio de grado N 1 que pasa exactamente por estos puntos. Este polinomio se puede expresar de la forma p(x) = c1 xN 1 + c 2 xN 2 + · · · + cN 1x + cN y verifica que p(xk ) = yk para k = 1,... , N , es decir: 8 > > y1 = c 1 xN 1 1 + c 2 xN 1 2 + · · · + cN 1 x1 + cN < y2 = c 1 xN 2 1 + c 2 xN 2 2 + · · · + cN 1 x2 + cN > >... : yN = c 1 xN N 1 + c 2 xN N 2 + · · · + cN 1 xN + cN La resolución de este sistema proporciona los valores de los coeficientes c1 , c2 ,... y cN. Este procedimiento se conoce como interpolación global de Lagrange. 1 y (x 2 , y 2 ) (x N , y N ) (x 1 , y 1 ) (x 3 , y 3 ) x Los valores de los coeficientes del polinomio se calculan habitualmente con ayuda de algún programa informático. En el ejemplo siguiente se explica cómo hacerlo con MATLAB. 1 Joseph Louis Lagrange (1736–1813), fue un matemático, físico y astrónomo italiano nacido en Turín, aunque vivió casi siempre en Francia y Rusia. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 178 Ejemplo 4.12 La temperatura del aire cerca de la tierra depende de la concentración K del ácido carbónico (H2 CO3 ) en él. En la tabla de más abajo se recoge, para diferentes latitudes L sobre la tierra y para el valor de K = 0.67, la variación K de la temperatura con respecto a una cierta temperatura de referencia. Calcular el polinomio de interpolación asociado a estos datos. L -11 -7 5 8 12 K -7 2 -3 4 -5 Aquí, la magnitud K es la variable dependiente, y L es la variable independiente: L !x K !y Se desean calcular, con MATLAB, los coeficientes del polinomio de grado 4 (ya que hay 5 datos) que toma dichos valores, es decir, encontrar un polinomio 8 > > p( 11) = 7 > > < p( 7) = 2 p(x) = c1 x4 + c2 x3 + c3 x2 + c4 x + c5 que verifique p(5) = 3 > > > > p(8) = 4 : p(12) = 5 Para ello, basta escribir las siguientes órdenes en MATLAB: x=[-11,-7,5,8,12] y=[-7,2,-3,4,-5] c=polyfit(x,y,4) Con esto se obtendrá: c = -0.0027 0.0048 0.3909 -0.2255 -10.5492 lo que significa que el polinomio interpolante es: p(x) = 0.0027 x4 + 0.0048 x3 + 0.3909 x2 0.2255 x 10.5492 La interpolación polinómica global no tiene mucho interés práctico (aunque sí lo tiene teórico), sobre todo cuando aumenta el número de datos que se quieren interpolar. Las razones principales son dos: Es inestable, es decir, una pequeña variación en los datos puede producir una gran diferencia en los polinomios de interpolación. Esto es muy importante cuando los datos proceden de mediciones, ya que es inevitable cometer errorres. Cuando aumenta el número de puntos a interpolar hay que recurrir a polinomios de grado cada vez mayor, y los polinomios de grados altos tienden a ser muy oscilantes, y normalmente no representan bien los valores de una función sin grandes variaciones. Este fenómeno se observa muy bien en el Ejemplo 4.13. Mucho más interés práctico tiene la interpolación a trozos, que se explica más adelante. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 179 Ejemplo 4.13 Para interpolar los valores: x = (0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15), y = (10, 20, 30, 10, 10, 10, 10.5, 15, 50, 60, 85) es necesario un polinomio de grado 10 (ya que hay 11 datos). Los puntos y el polinomio están representados en la figura siguiente: 150 100 50 0 −50 −100 −150 −200 −250 −300 −350 −400 0 2 4 6 8 10 12 14 16 Se observa que el procedimiento de interpolación global es, en general inestable, ya que los polinomios tienden a hacerse oscilantes al aumentar su grado y eso puede producir grandes desviaciones sobre los datos. Interpolación lineal a trozos Hablando en términos muy imprecisos, la interpolación lineal a trozos consiste en unir con segmentos rectos los pares de puntos consecutivos que se quieren interpolar. Consideramos N puntos (xk , yk ), k = 1,... , N , con los valores de xk todos diferentes y ordenados en orden creciente o decreciente. Se llama interpolante lineal a trozos a la poligonal que sobre cada intervalo formado por dos valores de x consecutivos [xk , xk+1 ], k = 1,... , N 1, está definida por el segmento que une los puntos (xk , yk ) y (xk+1 , yk+1 ), como en la Figura 4.20. 100 80 60 40 20 0 −20 0 2 4 6 8 10 12 14 16 Figura 4.20: Interpolante lineal a trozos. Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 180 Ejemplo 4.14 Con los mismos datos del Ejercicio 4.13, x = (0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15), y = (10, 20, 30, 10, 10, 10, 10.5, 15, 50, 60, 85) la interpolación lineal a trozos daría como resultado la función poligonal de la figura: 150 100 50 0 −50 −100 −150 −200 −250 −300 −350 −400 0 2 4 6 8 10 12 14 16 Compárese la diferencia de valores que se encontraría si se calculara el valor de la función y en x = 1 con cada uno de los interpolantes: con el interpolante polinómico del Ejercicio 4.13 se obtendría el valor y = 247.0336, mientras que el interpolante lineal a trozos se obtendría y = 10. Ajuste de datos La técnica de interpolación que hemos explicado antes requiere que la función que interpola los datos pase exactamente por los mismos. En ocasiones esto no da resultados muy satisfactorios, por ejemplo si se trata de muchos datos. También sucede con frecuencia que los datos vienen afectados de algún error, por ejemplo porque provienen de mediciones. No tiene mucho sentido, pues, obligar a la función que se quiere construir a «pasar» por unos puntos que ya de por sí no son exactos. Otro enfoque diferente es construir una función que no toma exactamente los valores dados, sino que «se les parece» lo más posible, por ejemplo minimizando el error, medido éste de alguna manera. Cuando lo que se minimiza es la suma de las distancias de los puntos a la curva hablamos de ajuste por mínimos cuadrados. La descripción detallada de este método se escapa de los objetivos de estas notas. En el siguiente Ejemplo se muestra cómo calcular con MATLAB la recta y la parábola que mejor se ajustan a unos datos. (0.9, 0.9) (1.5, 1.5) (3, 2.5) (4, 5.1) (6, 4.5) (8, 4.9) (9.5, 6.3) Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla 4. Métodos numéricos 181 Ejemplo 4.15 Se desea calcular la recta y la parábola que mejor se ajustan a los datos siguientes: x 0.9 1.5 3 4 6 8 9.5 y 0.9 1.5 2.5 5.1 4.5 4.9 6.3 1. Cálculo de la recta y = ax + b que mejor se ajusta a los siguientes datos. Dicha recta se llama recta de regresión. En MATLAB, escribir las órdenes siguientes: x=[0.9 , 1.5 , 3 , 4 , 6 , 8 , 9.5] y=[ 0.9 , 1.5 , 2.5 , 5.1 , 4.5 , 4.9 , 6.3] c=polyfit(x,y,1) Se obtendrá c = 0.5688 0.9982 ⇡ (0.57, 1) lo que significa que la recta buscada es y = 0.57 x + 1 2. Si lo que se desea es calcular la parábola de regresión: x=[0.9 , 1.5 , 3 , 4 , 6 , 8 , 9.5] y=[ 0.9 , 1.5 , 2.5 , 5.1 , 4.5 , 4.9 , 6.3] c=polyfit(x,y,2) Se obtendrá c = -0.0617 1.2030 -0.0580 ⇡ ( 0.06, 1.2, 0.06) lo que significa que la parábola buscada es y = 0.06 x2 + 1.2 x 0.06 7 6 5 4 3 2 1 0 −1 0 1 2 3 4 5 6 7 8 9 10 Matemáticas Aplicadas a la Biología - Grado en Biología R. Echevarría - Dpto. EDAN - Univ. de Sevilla