Source: https://www.scribd.com/document/55907017/mastquim
Timestamp: 2018-03-17 17:06:33+00:00

Document:
Uploaded by Cesar David Florez
Métodos numéricos para las ecuaciones diferenciales Aplicaciones a la Química
Jose S. Cánovas Peña 5 de febrero de 2010
Índice General Advertencia: Esta es la primera versión de los apuntes de métodos numéricos del cuarto curso de Ingeniería Industrial. No han sido corregidos y probablemente contengan numerosos errores, pero he decidido colgarlos en la web para que sirva de guia de estudio. Para corregir errores, podeis escribir a mi dirección de correo Jose.Canovas@upct.es.
.5 Extrapolación de Richardson . . . . . . . . 1. . . . . .1 Introducción . . . . . . . . . . . . . . .2 El método de cuatro etapas . . . . . . . . . . . . . . . . . . . 2. . 1 3 3 4 6 8 8 9 12 13 14 14 16 17 18 23 23 23 25 27 29 33 35 35 36 37 38 39 40 43 45 1 Introducción a las ecuaciones diferenciales: modelos de la Química 1.1 Error local de truncamiento en el método de Taylor . . . . . . . . . . 1. . . . .2 Error local de truncamiento en los métodos de Runge—Kutta . . . . . . . . . . . 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Métodos de un paso 2. . . . . . . . . .1 El método de 3 etapas . . . . . . . . . . . . . . 1. . .4. . . . . . . . . . . .2 Climatización de ediﬁcios con varias estancias . . . . . . . . .1 Ecuaciones diferenciales . . . . . .4 Relación entre el error local y el error local de truncamiento . . . . . . .4. . . . . . . . . . . . . . . 2.2 Ley de enfriamiento de Newton. . . . . . . . . . . . . . . . . . . .6. . . . . . . . . . . . . . . . . .4. . . . . . . 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. . . . . . . . . . . 2. . . . . . . . . . . . . . . . . . 2. . . . . . . . . . . . . . 2. . . . . . . . . . . . .3 Métodos de Runge—Kutta de orden dos . .4. .4 Análisis del error en los métodos de orden uno . . . .2. . . . . . . . . . . . . . . . . . .5. . . . . . .4. 1. . . . . . . . . . . .3 Sistemas de ecuaciones diferenciales . . . . . . .2 Método de Taylor de orden 2 . . . . . . . . . . . . 1. . . .6 Más sobre los métodos Runge—Kutta .2 Métodos de Taylor . . 1. . . . 1. . 1. . . . . . . . . . . .2 Soluciones de ecuaciones diferenciales . . . . . 2. . .4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3 Métodos de más etapas . . . . . . . . . . . . . . . . . . . . . . . . 1. . . . . 2. . . . . . .6. . . . .3 Cinética de las reacciones químicas . . . . . . . . . . . . . . . . . . . . . . . . .4 Modelos de la química descritos por una ecuación . . .4 Cinética de las reacciones químicas . . . . . . . . . . . . . . .5. . . . . . . . . . . . . . . . . . . . . 2. . . . . . . . . . . . . . . .1 Método de Euler . . . . . . . . . 2. . . . .3 Relación entre el error local de truncamiento y el error global . . . . . . . . 2. . . . . . . . . . .5 Modelos de la química descritos por un sistema de ecuaciones . . . . . . . . . . . . . . . . . . 2. . . . . . . . . . . . . . . .5. . . . . . . . . . . . . . . .1 Problemas de mezclas con varios recipientes . . .4. . . . . . . . .1 Descomposición radioactiva . . . . . . . . .4. . . . . . . .6.Índice General I Bloque de teoría . 1. . . .2. . . . . . . . iii . . . . . . . . 1. . . . . . . . . . . . . . . . . . . . 2.6 Ejercicios .3 Problemas de mezclas químicas. . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . 3. . . . . . . . . . . . . . . . . .2. . . . . . . . . . . 4. . . . . . . . . . . . . . . . . . . . . .2. . .4 Do . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . constantes.2 Operaciones . . . . . . . . . . . . . . 6 Programación en Mathematica 6. .7 Metodos predictor—corrector . . 4. . . . . . . . . . .2. . .1 Operaciones y deﬁnición de variables. .3. . . . . . .3 Aplicaciones de las ecuaciones lineales con coeﬁcientes 5. . . 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2 Which . . . . . . . . . . . . . .2 Movimiento amortiguado. . . . . . 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3 Movimiento forzado. . .1 Deﬁnición de variables y funciones . . . .4 Aplicación a los circuitos eléctricos . . . . . . . .3 Funciones matemáticas de interés . . . . . 6. . . . . . . . . . . . . . . . . . . . . 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. . . . . . . . . . . . . . . .5 Estabilidad de los métodos multipaso lineales . . . . . . . . . . . . . . . . .4 Métodos de multipaso deducidos a partir de la integración numérica 3. . . .1 Introducción . . . . . . . . . . . . . . corchetes y llaves . . . . . . . . . . . . . . . .2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. . . . . . . . . . . . 6. .1 Preliminares . . . . . . .4 Aprovechando cálculos anteriores . . . . . . . . .2 Paréntesis. . . . . . . . . . . . . . 3. . . . . . . . .2. . . .2 Ecuaciones diferenciales lineales. . . . . . . . . . . . 5. . 5. . . . . . . . . .3. . . . .1. 6. . . . . . .1 Movimiento armónico simple. . . . . . . . . . . . . . . . . . . . 6. . . . . . . . . . . . . .3 Primeros ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 63 63 65 66 66 67 68 70 71 73 73 75 76 76 77 78 78 81 81 81 83 84 84 85 85 86 86 87 87 4 Introducción a Mathematica 4. . . . . . . . . .2. . . .7 Break . . . . . . . . . . . . .3. . . . . . . .6 For . . . . . . . . II Bloque de práctica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. . . . . . . . .Índice General 3 Métodos multipaso 3. . . . . . . . . . . . .6 Derivadas de funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4. . 4. . . . . . .1 Métodos de Adams—Bashforth . . . . . . . . . .2. . . . . . . . . . . . . . 3. . . . 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 If . . . . . . . 6. . . . . . . . . . . . . . 6. . . . . . . . . . . . . . . . . . . . . . . 5. . .1 Ecuaciones diferenciales de primer orden . . . . . . . . . . . . . . . . . 5. . . . . . .7 Representación gráﬁca de funciones 5 Ecuaciones diferenciales con Mathematica 5. . 47 47 47 49 51 52 53 55 57 58 59 . . . . . . . . . . . .1 Errores . . .5 While . . . . . . . . . 6. 6. . .1. . 6. .3 Switch . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2 Construcción de programas con Mathematica. . . . . . 6. . . 4.4. .8 Multipaso o Runge—Kutta . . . . . 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .6 Fórmulas BDF . . . . . . . . . . . . . . .5 Deﬁnición de variables y funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2 Métodos de Adams—Moulton . . . . . . . iv . . . . . . . . . . . .2 Métodos multipaso lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2. . . . . . . . . . . . .
4 Fórmula de inversión compleja . . .2. . . . . . .2. . . . .2 Transformada Z inversa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2. A. . . . .2. . . . . v . . . . . . . . . . . . . . . . . . . . A. . 87 88 88 89 90 95 95 98 98 100 101 103 105 106 A Ecuaciones en diferencias A. . . . . . . . . . . . . . . . . . . . . . .3 Ecuaciones no lineales . . . . . . . . . . . . . . . . . . . . . . . . . .3 Aplicación a la resolución de la ecuación en diferencias lineales con coeﬁcientes constantes . . . . . . . . . . . . . . . . . . . 6. . . .10 Input . . . . . . . . . . . . . . . . . . .2. . . . . . .8 Goto. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3 Un programa . . . . . . . . . . . . . . 6. . . . . . . . . . . . . . .4 Presentaciones gráﬁcas . . . . . . . . . . . . . . .Índice General 6. . . . . . . A.2. . . . . . .2. . . . . . . . . . . . . . . . . . . . . 6. . . . .1 Introducción . . . . . . . . A. . . .2 Transformada Z . . . A. . . . . .5 Funciones de transferencia. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1 Deﬁnición y propiedades básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. . . . . . . 6. . . . . . . .9 Print . . . . A. . . . . . . . . . . . . . . . . . . Label . . . . . . . . . . . . . . . . . . . . . . . . .
Parte I Bloque de teoría 1 .
b). Problemas de mezclas. que tienen órdenes 2. b). Ley del calor de Newton. y(x).1) al valor de la derivada más alta en dicha expresión. . Concepto de solución. y(x). Como vemos. . y 0 .1) si existe la derivada n—ésima de y en todo punto del intervalo (a. 3 (1. y. que llamaremos variable independiente y las n primeras derivadas respecto de x de una variable y. (x. Ecuación y sistemas de ecuaciones diferenciales.1) Una ecuación diferencial una expresión de la forma donde F es una función real deﬁnida en un cierto abierto A ⊆ Rn+2 . Descomposición radioactiva.2) . 1.. que se llama variable dependiente por ser una función dependiente de la variable x.. y · y 0 · y 00 = x.. y 0 (x). Diremos que una función y : (a.1 Ecuaciones diferenciales F (x. . una ecuación diferencial es una expresión en la que aparecen ligadas una variable x.. Ejemplos de ecuaciones diferenciales son los siguientes: y 00 + log(xy) − x = y. b) y F (x. Cinética de las reaciones químicas. (1. y 3) + xy 0 + ex sinh y = 0. y n) ) = 0.... 3 y 2. y n) (x)) ∈ A para todo x ∈ (a. y n) (x)) = 0 para todo x ∈ (a. e y(x) es una función real de variable real. Por ejemplo tomemos la ecuación diferencial de orden uno y 0 − y tan x = 0. b) ⊆ R → R es solución de la ecuación (1. Se llama orden de la ecuación (1. y 0 (x)..Capítulo 1 Introducción a las ecuaciones diferenciales: modelos de la Química Sumario. respectivamente..
. c c sin x − tan x = 0. dada una ecuación diferencial de la forma F (x. de la forma F (y.. a0 y b funciones reales de variable real. cn ) = 0 (1. donde f : A ⊆ Rn+1 → R. y. En efecto. y n) ) = 0..2) las soluciones de una ecuación diferencial en caso de existir no son únicas.4) (1..3) con ci ∈ R. y las ecuaciones lineales a0 (x)y n) + a1 (x)y n−1) + . donde F no depende de la variable independiente explícitamente aunque ésta se encuentre impícita en la función y. . las soluciones de la misma pueden escribirse como g(x.... 2. a1 . y. Estas ecuaciones son obtenidas cuando sea posible despejar y n) en (1. Así las soluciones de una ecuación diferencial de orden n generan una familia n—paramétrica de curvas en el plano. En el ejemplo anterior. y 0 (x)) ∈ A para todo x ∈ (−π/2. c1 . y(x). esta función es una vez derivable con derivada c sin y 0 (x) = cos2 x . como puede verse fácilmente derivando de forma implícita la expresión x2 + y 2 = c respecto a la variable independiente x.. A lo largo del curso estudiaremos fundamentalmente ecuaciones resueltas respecto de la derivada de mayor orden de la ecuación.. an−1 . y. El dominio de deﬁnición de F es en este caso A = {(−π/2 + 2kπ.1). y además satisface que x y 0 (x) − y(x) tan x = para todo punto x ∈ (−π/2.. Así las curvas x√+ √ = c > 0 0 deﬁnen implícitamente las soluciones de la ecuación yy + x = 0 deﬁnidas en (− c. . y. donde c es una constante arbitraria es una solución de dicha ecuación diferencial. se veriﬁca que (x.. y n) ) = 0. Entonces c la función y : (−π/2.. π/2).. + an−1 (x)y 0 + an (x)y = b(x) con an .. es decir. . Con frecuencia las soluciones de la ecuación (1. y n−1) ).1) no podrán obtenerse de forma explícita y vendrán 2 y2 dadas de forma ímplicita por una ecuación de la forma g(x. y 0 . π/2 + 2kπ) : k ∈ Z} × R2 .. 2x cos cos x 1. sino que dependen de ciertas constantes arbitrarias provenientes de la integración.. i = 1. y 0 ... y) = 0. . π/2) → R dada por y(x) = cos x .. .2 Soluciones de ecuaciones diferenciales Como vimos en el ejemplo (1. la solución y(x) = c/ cos x deﬁne 4 . ecuaciones de la forma y n) = f (x. n. π/2). y 0 ) = y 0 − y tan x. En general. . y 0 . Serán también de especial interés para nosotros las ecuaciones autónomas..Introducción a las ecuaciones diferenciales Esta ecuación viene deﬁnida por la función F : A ⊂ R3 → R dada por F (x. c).
c1 .4) puede construirse una ecuación diferencial de la manera siguiente.. Si consideramos el problema de condiciones iniciales ½ 0 y − y tan x = 0 y(0) = 1. n−1) Lo que se espera añadiendo estas condiciones es eliminar los parámetros de la familia n—paramétrica de soluciones. y0 .. ⎪ ⎪ ⎪ y(x0 ) = y0 . obtendremos una ecuación diferencial de orden n dada por (1. ... que no depende de parámetro alguno. eliminando los parámetros c1 . . Derivando n veces (1. Mención aparte merecen aquellas que no aparecen comprendidas en la familia n—paramétrica. Despejando c1 y c2 y sustituyendo en la primera ecuación tenemos que y 00 = y es la ecuación diferencial que deﬁne a la familia de curva anteriores. la ecuación y 2 + (y 0 )2 = −1 no posee ninguna solución. y 0 .. y 0 . Por ejemplo.. 5 . y n) ) = 0 0 donde las consde n condiciones iniciales de la forma y(x0 ) = y0 . . Por ejemplo y 0 = −2y 3/2 tiene como solución general y(x) = 1/(t + c)2 y como solución singular y(x) = 0. Además la ecuación de orden uno (y 0 − y)(y 0 − 2y) = 0 tiene por soluciones a las funciones dadas por la expresión (y − c1 et )(y − c2 e2t ) = 0. a partir de una familia n—paramétrica de curvas deﬁnida por (1. ⎪ . las soluciones de las ecuaciones diferenciales vienen dadas por una familia n—paramétrica de curvas y por lo tanto la solución no es en general única.. derivando implícitamente respecto de x tenemos que y 0 = c1 ex − c2 e−x . de manera que encontremos la solución de la ecuación diferencial satisfaga adicionalmente estas condiciones. Nótese que es una ecuación de orden dos dado que la familia depende de dos parámetros... Se deﬁne un problema de condiciones iniciales o de Cauchy al problema de la forma ⎧ ⎪ F (x. π/2) → R dada por y(x) = cos x .y n−1) (x0 ) = y0 n−1) 0 tantes x0 . . c2 ∈ R. que c recordemos. Nótese que esta deﬁnición es ambigua y depende de la familia de curvas al ser y(x) = C 2 /(Cx + 1)2 una familia uniparamétrica de soluciones de y 0 = −2y 3/2 conteniendo la solución nula. Nótese que se añaden tantas condiciones iniciales como orden tiene la ecuación. Por ejemplo. mientras que y 2 + (y 0 )2 = 0 tiene como única solución y(x) = 0. Recíprocamente.. las llamadas soluciones singulares. y. De acuerdo con lo visto anteriormente. y0 son números reales..3).2).. Las soluciones obtenidas como familia n—paramétrica de curvas se llaman soluciones generales de la ecuación diferencial. Para evitar este hecho. tomemos la ecuación de orden uno (1. Sin embargo...4) respecto de x obtenemos n + 1 ecuaciones de las que. Por ejemplo. ⎪ ⎨ 0 0 y (x0 ) = y0 .. obteniendo entonces una solución que sea única. y n) ) = 0. y 00 = c1 ex + c2 e−x . suele acompañarse a una ecuación diferencial F (x. si consideramos la familia de las curvas del plano dependiente de dos parámetros dada por la ecuación y = c1 ex +c2 e−x . ⎪ . y. y 0 (x0 ) = y0 . ⎪ .Introducción a las ecuaciones diferenciales una familia de curvas en el plano dependiente del valor o parámetro de c. y0 . ⎪ ⎪ n−1) ⎩ n−1) y (x0 ) = y0 .. c2 . cn . tenía por solución y : (−π/2. donde c es una constante arbitraria. no siempre las ecuaciones diferenciales de orden n presentan una solución general que se expresa mediante familias n—paramétricas.
⎪ ..Introducción a las ecuaciones diferenciales tendríamos que necesariamente 1 = y(0) = c/ cos(0) = c.. ⎪ . y 0 . y1 .. y2 . y .. Sólo para el caso de los sistemas de ecuaciones diferenciales lineales con coeﬁcientes constantes. y 0 . ⎩ 2 0 y3 = y1 y2 y3 . y1 . 0 y2 = x + y1 + y2 . Primero claro está. Estamos interesados en aquellos sistemas de ecuaciones diferenciales en los que podemos despejar la primera derivada de cada una de las funciones incógnita. y2 ... ym ).3 Sistemas de ecuaciones diferenciales donde y1 . ⎪ 0 ⎩ y = f (x. 1 2 m . que veremos un poco más tarde existen algoritmos que permiten el cálculo explícito de las soluciones. ⎪ ⎪ 0 ⎨ y = f2 (x. ⎧ 0 2 ⎨ y1 = xy1 + y2 − y3 . . y 0 ) = 0. salvo en casos excepcionales. . 1 ≤ i ≤ m.... y 0 . ym ). .. y2 . . debemos deﬁnir qué 6 donde fi : A ⊆ R1+m → R. 1 ≤ i ≤ m. y1 .. ⎪ .. 1. y2 . ⎪ ⎩ F (x. y . ninguna puede deducirse de las demás. . por lo que c = 1 y la única solución del problema de condiciones iniciales es y(x) = 1/ cos x. y2 . y2 . .. .. m 1 1 2 2 m m En general la resolución de estos sistemas no es posible.. y . m m 1 2 m Un sistema de ecuaciones diferenciales es una expresión de la forma ⎧ 0 0 0 ⎪ F1 (x.. y1 . Sin ir más lejos. y 0 . Se suele suponer que hay igual número de ecuaciones que de incógnitas de manera que todas las ecuaciones son independientes. 0 y = x + y1 + y2 y3 . y ). Sin embargo.. y1 . . Ejemplos de estos sistemas son ½ 0 2 y1 = xy1 + y2 . es relativamente sencillo saber cuándo un sistema tiene solución. son funciones reales de varias variables. y . sistemas de la forma ⎧ 0 ⎪ y1 = f1 (x. son funciones reales. ym son funciones reales a determinar que dependen de x y Fi : A ⊆ R1+2m → R. es decir. ym ) = 0. ym . ⎪ ⎪ ⎨ F2 (x. y 0 ) = 0. es decir.. ym .. 2 .. o más precisamente cuándo un problema de condiciones iniciales asociado tiene solucón. ⎪ . el problema ½ 2 y + (y 0 )2 = 0 y(0) = 0 no tiene solución y ½ y 0 = 3y 2/3 y(0) = 0 tiene al menos dos soluciones dadas por y(x) = 0 e y(x) = x3 . y .. Sin embargo esta estrategia no siempre produce los frutos deseados.
. 0 ⎪ y3 = y1 y2 y3 . y . y1 . ym ).. y3 (0) = 1. son funciones reales continuas en el ∂fi abierto A. y3 (0) = 1. es decir. .0 ⎪ y = f (x. ⎪ 0 ⎨ y2 = x + y1 + y2 y3 . y .. . y1 . es un problema de condiciones iniciales.. . . y1 . Theorem 1 Sea el problema de condiciones iniciales ⎧ 0 ⎪ y1 = f1 (x. ⎩ y1 (0) = 2. y2 (0) = 0. y . y1 .. y1 . ⎪ ⎪ 0 ⎪ y = f2 (x.. donde x0 . y2 .. ⎪ 0 ⎪ y3 = y1 y2 y3 .. ⎪ 2 ⎨ . 7 . Por ejemplo ⎧ 0 2 ⎪ y1 = xy1 + y2 − y3 . y (x ) = y 1 0 1 2 0 2 m 0 m donde (x0 .. y ). y (x ) = y 1 0 1 2 0 2 m 0 m junto con las condiciones yi (x0 ) = yi .. aunque pensamos que ésta es de un nivel excesivo para alumnos de primer curso.. cuya prueba puede verse en [Jim]. no sería un problema de condiciones iniciales. . ym ) ∈ A... Por ejemplo el problema que consideramos anteriormente ⎧ 0 2 ⎪ y1 = xy1 + y2 − y3 . Nótese que todas las condiciones iniciales implican el conocimiento de la función en 0. ⎪ ⎩ y1 (0) = 2. Supongamos además que las funciones ∂yj existen y son continuas en A.. ⎪ ⎩ y1 (0) = 2. . ..Introducción a las ecuaciones diferenciales entendemos por un problema de condiciones iniciales para sistemas de ecuaciones diferenciales. Para el caso de los problemas de condiciones iniciales para sistemas de ecuaciones diferenciales tenemos el siguiente resultado análogo al de ecuaciones diferenciales de orden uno.. y (x ) = y ... y (x ) = y .. ym son números reales. y2 ... ym ). . fi : A ⊆ R1+m → R. deﬁnido en un intervalo abierto I de la recta real. y ). ⎪ m ⎪ m 1 2 m ⎪ ⎩ y (x ) = y . ⎪ ⎪ 0 ⎪ y = f2 (x. y2 . ym ). 1 ≤ i ≤ m. .. y . y2 . ⎪ ... Entonces existe una solución del problema de condiciones iniciales anterior yi : I → R. y3 (0) = 1.. ya que conocemos y2 en 1 e y1 e y3 en 0. . ym ). ⎪ 0 ⎨ y2 = x + y1 + y2 y3 . 0 ⎪ y3 = y1 y2 y3 .. ⎪ 0 ⎨ y2 = x + y1 + y2 y3 . y1 . 1 ≤ i ≤ m. y2 (1) = 0. lo siguiente ⎧ 0 2 ⎪ y1 = xy1 + y2 − y3 ... ⎪ . y2 (0) = 0. Este resultado es fácil de aplicar. ⎪ m ⎪ m 1 2 m ⎪ ⎩ y (x ) = y . ⎪ 2 ⎨ .0 ⎪ y = f (x.. y2 . . .. Dicho problema es un sistema de ecuaciones diferenciales ⎧ 0 ⎪ y1 = f1 (x.
. ⎪ ⎪ ⎪ y(x0 ) = y0 . ⎪ n 1 2 n ⎪ ⎪ ⎩ n−1) 0 y1 (x0 ) = y0 . y1 . 8 . es decir. continuas y las derivadas parciales de cada función respecto de y1 . y2 = y 0 .4 1. y ). y2 (x0 ) = y0 . Una sustancia radioactiva tiende a descomponerse con el tiempo formando nuevas sustancias y liberando a su vez una gran cantidad de energía. y1 . . satisface la ecuación diferencial y 0 = Ky donde K es una constante que depende de la sustancia considerada. es decir. yn (x0 ) = y0 .1 Modelos de la química descritos por una ecuación Descomposición radioactiva Consideremos un isótopo radioactivo del cual tenemos una cantidad y(t) que varía con el tiempo t. yn . ⎪ ⎨ 0 0 y (x0 ) = y0 .. . 0 ⎪ yn−1 = yn .. y3 = y 00 . ⎪ . con las que el sistema anterior queda como ⎧ 0 ⎪ y1 = y2 . .Introducción a las ecuaciones diferenciales 2 es tal que f1 (x. f2 (x. Entonces este problema de condiciones iniciales tiene solución única. y .. y proporciona las soluciones y(t) = CeKt . y2 .4.. y2 e y3 son continuas. Esta ecuación es de variables separadas. con C la constante proveniente de la integración. El Teorema 1 tiene una aplicación inmediata a los problemas de condiciones iniciales de orden n. ⎪ ⎪ 0 ⎪ y = f (x.. Por ejemplo. la función y y sus n − 1 primeras derivadas. y introduzcamos las variables y1 = y. y1 . ⎪ ⎪ ⎪ . ⎪ 0 ⎪ ⎪ y2 = y3 .. .. y 0 . y2 . Se ha comprobado experimentalmente que la velocidad con que una sustancia radioactiva se descompone es directamente proporcional a la cantidad de sustancia existente en dicho instante. ⎨ . y3 ) = y1 y2 y3 son funciones deﬁnidas en R4 .. y . aunque no tengamos ni idea de cómo calcularla. ⎪ . y3 ) = x + y1 + y2 y3 y f3 (x.... ...yn = y n−1) . y3 ) = xy1 + y2 − y3 . y2 . y. 1. ⎪ . y y entonces el Teorema 1 nos garantiza la existencia y unicidad de soluciones siempre que la función f sea continua y de clase C 1 respecto de las variables y1 . ⎪ ⎪ n−1) ⎩ n−1) (x0 ) = y0 . consideremos el problema ⎧ n) ⎪ y = f (x. y n−1) )..
como el tiempo necesario para que una cantidad de dicha sustancia se reduzca a la mitad.Introducción a las ecuaciones diferenciales Se deﬁne la vida media de una sustancia radioactiva tm . K Como podemos observar. Obtenemos la solución de la ecuación diferencial. ¿cuándo podremos tomar el café si la temperatura idónea para tomarlo es de 650 C ? Para responder a esta pregunta. 9 . 20 − y(t) que nos proporciona la solución y(t) = Ce−Kt + 20. si servimos una taza de café a una temperatura de 950 C y al minuto está a 850 C. la vida media de la sustancia depende de la constante K. donde K es una constante proveniente de la integración. su vida media puede calcularse resolviendo primero el problema de condiciones iniciales ¾ y 0 = Ky y(0) = N0 que proporciona la solución y(t) = N0 eKt . ésta veriﬁca la ecuación diferencial y 0 = K(T − y). Si tenemos una cantidad inicial de una sustancia N0 . Supongamos que tenemos un cuerpo inerte que no produce calor de manera autosuﬁciente.2 Ley de enfriamiento de Newton. basta con resolver el problema de condiciones iniciales ¾ y 0 = K(20 − y) y(0) = 95. De observaciones experimentales se sabe que la temperatura superﬁcial de dicho cuerpo varía con una rapidez proporcional a la diferencia de temperatura entre el objeto y su entorno. calculamos dicha constante resolviendo la ecuación 95 = y(0) = C + 20.4. y suponiendo que la habitación está a 200 C. si denotamos por y(t) la temperatura del cuerpo con el tiempo. Por ejemplo. que es intrínseca de cada sustancia. cómo por ejemplo el agua. y posteriormente la ecuación N0 = N0 eKtm . una piedra o un reloj. que es de variables separadas calculando Z Z y 0 (t) dt = Kdt. tm = − 1. Es decir. Los contenidos de esta sección pueden verse en [NaSa]. Imponiendo ahora que y(0) = 95. 2 con lo que la vida media es log 2 . donde K es una constante de proporcionalidad y T es la temperatura ambiente en ese momento. Dicho comportamiento es conocido cómo la ley de enfriamiento de Newton.
Si el ediﬁcio no esta vacío se produce un calentamiento adicional debido al calor corporal. Un valor normal para un ediﬁcio cerrado oscila entre la 2 y las 4 horas para la constante 1/K. log(13/15) Es decir. Si suponemos constante E(t) = E0 . como al minuto de haber servido el café la temperatura de éste había descendido hasta los 850 C tenemos que 85 = y(1) = 75e−K + 20. que escribiéndola como vemos claramente que es lineal. máquinas en funcionamiento.57 minutos. aproximadamente unos tres minutos y medio después de haber servido el café. El tiempo que transcurre desde el valor T (0) − E0 hasta el valor (T (0) − E0 )/e es t0 = 1/K. la ecuación anterior queda como T 0 (t) = K(E(t) − T (t)) + H(t) + U(t).Introducción a las ecuaciones diferenciales con lo que C = 75. cuya razón denotaremos por H(t). luces. Si adicionalmente el ediﬁcio dispone de un sistema de calefacción o de aire acondicionado. que recibe el nombre de constante de tiempo del ediﬁcio. Por otra parte. t= Aplicación a la climatización de ediﬁcios Supongamos que tenemos un ediﬁcio que en un principio vamos a considerar como una unidad. d dT (t) = (T (t) − E0 ) = −K(T (t) − E0 ). Para averiguar el momento en el cual la temperatura de dicha taza es de 650 C basta resolver la ecuación 65 = 75et log(13/15) + 20. T 0 (t) = −KT (t) + (KE(t) + H(t) + U(t)) 10 T (t) − E0 = ce−Kt . que permite obtener el valor de la constante K = log(13/15). dt dt . es decir. etcétera. entonces la ecuación puede escribirse como T 0 (t) = que nos proporciona la solución Si T (0) es la temperatura inicial del ediﬁcio c = T (0) − E0 y la solución es T (t) − E0 = (T (0) − E0 )e−Kt . Entonces. Si T (t) es la temperatura del ediﬁcio vacío en un instante de tiempo t y E(t) es la temperatura en el exterior (que puede ser variable). la ley de Newton aﬁrma que T 0 (t) = K(E(t) − T (t)). se produce un aumento o disminución de la temperatura que denotaremos por U(t). c ∈ R. no vamos a tener en cuenta el número de habitaciones que tiene (ya veremos posteriormente este caso). y que suele medirse en horas. La función y(t) = 75et log(13/15) + 20 deﬁne entonces la evolución de la temperatura de la taza de café con el tiempo. que da la solución log(9/15) = 3.
9 o C. A mediodia se apaga el aparato de aire acondicionado y la gente se va a sus casas. 11 T (t) = −246e−t/64 + 276. 27 = −15e−t0 /4 + 35. Si la constante de tiempo del ediﬁcio es de 4 horas. ¿cuál será la temperatura del ediﬁcio a las 2 de la tarde? ¿En que momento la temperatura en el interior será de 27 o C? Para responder a esta pregunta planteamos la ecuación diferencial 1 T 0 (t) = (35 − T (t)). El tanque se encuentra bien aislado y tiene una constante de tiempo de 64 horas.Introducción a las ecuaciones diferenciales Example 1 Supongamos una mañana de sábado caluroso que en una tienda. La temperatura exterior permanece constante a 35 o C. La solución de dicha ecuación diferencial es T (t) = ce−t/64 + 276. 4 dado que H(t) = U (t) = 0. junto con la condición inicial T (0) = 20. . El momento t0 en que la temperatura será de 27 o C se obtendrá al resolver la ecuación que nos da 8 ' 2. t0 = −4 log T (t) = −15e−t/4 + 35. 64 junto con la condición inicial T (0) = 30 o C. que se corresponde con la temperatura al mediodía.51 horas 15 es decir. El panel solar genera 2000 kilocalorías por hora durante el día y el tanque tiene una capacidad caloríﬁca de 2 o C por cada 1000 kilocalorías. La solución de la ecuación diferencial será T (t) = ce−t/4 + 35. aproximadamente a la 2 horas y media. de donde la solución del problema de condiciones iniciales es Al cabo de 12 horas la temperatura del agua del tanque es T (12) = 72. mientras las personas están trabajando el aire acondicionado mantiene la temperatura de la tienda a 20 o C. U(t) = 2 o C/1000 Kcal × 2000 Kcal/h = 4 o C/h. Si el agua se encuentra inicialmente a 30 o C y la temperatura ambiente es de 20 o C. y con la condición inicial obtenemos c que nos proporciona la solución Así a las dos de la tarde la temperatura será de √ T (2) = −15/ e + 35 ' 25. ¿cuál será la temperatura del tanque al cabo de 12 horas de luz solar? En este caso con lo que la ecuación diferencial que modeliza el fenómeno es 1 T 0 (t) = (20 − T (t)) + 4.06 o C. Example 2 Un calentador solar de agua consta de un tanque de agua y un panel solar.
Vamos a determinar la cantidad de sal que habrá en el recipiente al cabo de 20 minutos. Por el método de variación de constantes calculamos la solución de la ecuación no homogénea imponiendo que y(t) = K(t)et/400 sea solución de la misma.4 Z et/400 et/400 = 0. donde ve y vs son las velocidades de entrada y salida de dicha sustancia respectivamente. 12 . tenemos que la variación de dicha cantidad viene dada por y 0 = ve − vs . teniendo en cuenta que el volumen se mantiene constante. el problema de condiciones iniciales ¾ y y 0 = be − V0 +et−f t f y(0) = X0 modeliza la cantidad de sustancia que hay en el recipiente por unidad de tiempo. 4 + K(t) . Además sacamos parte de la solución que se produce dentro del recipiente a una velocidad de f m3 /sg. y que introducimos otra solución que contiene una concentración b Kg/m3 de dicha sustancia la cual es introducida en el recipiente a una velocidad de e m3 /sg. supongamos una tanque que contiene originalmente 400 litros de agua limpia. y se deja que la mezcla salga del recipiente a la misma rapidez. Entonces K 0 (t)et/400 + K(t) con lo que K(t) = 0. Supongamos que inicialmente teníamos una cantidad de X0 kilogramos de una sustancia diluida en una concentración de X0 /V0 Kg/m3 . La ecuación homogénea y 0 = 400 tiene por solución y(t) = Ket/400 . como y(0) = 0. En estos problemas aparecen involucradas sustancias. Como ve = be Kg/sg y vs = f · y(t)/V (t) Kg/sg donde V (t) = V0 + et − f t es el volumen de disolución en el recipiente por unidad de tiempo.3 Problemas de mezclas químicas.Introducción a las ecuaciones diferenciales 1.05 kilogramos de sal por litro a una velocidad de 8 litros por minuto. y La ecuación diferencial implicada es lineal.4. Además. 400 400 e−t/400 dt = −160e−t/400 + C. donde K es la constante procedente de la integración. Las ecuaciones diferenciales también tienen aplicación dentro de los problemas de mezclas. Así la solución de la ecuación diferencial será y(t) = −160 + Cet/400 . las cuales se mezclan dentro de un recipiente de volumen dado V0 . Si denotamos por y(t) la cantidad de sustancia en cuestión dentro del recipiente por unidad de tiempo. Para ello. planteamos el problema de condiciones iniciales ¾ y y 0 = 0. tenemos que 0 = −160 + C.4 − 400 y(0) = 0. Vertemos en el tanque agua que contiene 0. Por ejemplo.
y la solución del problema de condiciones iniciales es y(t) = 160(et/400 − 1). La ley de velocidad diferencial establece que − 1 d[A] 1 d[C] = = k[A]n [B]m . a dt c dt será donde n y m son enteros positivos o la mitad de números enteros positivos. supongamos que incialmente hay cuatro y dos moles/litro de cada reactivo. Vamos a determinar la cantidad de producto C que tendremos a las 2 horas. No obstante. a dt b dt c dt d dt d[C] dt Nótese que d[A] será negativo al estar desapareciendo el reactivo mientras que por ejemplo dt positivo. y que al cabo de uno hora tenemos un mol/litro de C. La constante k recibe el nombre de constante de reacción. Si denotamos por [A] la concentración en moles por litro del reactivo A. la cantidad de sal que hay dentro del tanque es y(20) = 160(e1/20 − 1) ' 8. dt dt dt 13 . es decir. como NO + O3 → NO2 + O2 . hay reacciones que se llaman elementales. Denotemos por x(t) la concentración de C en cada instante.Introducción a las ecuaciones diferenciales con lo que C = 160. Dado que − d[A] d[B] d[C] =− = . y en las que los número n y m son iguales a uno. se veriﬁca que las velocidades de descomposición o formación de cada uno de los elementos de la reacción satisfacen la relación − 1 d[B] 1 d[C] 1 d[D] 1 d[A] =− = = . esquematizado como aA + bB → cC + dD. A los 20 minutos.4.4 Cinética de las reacciones químicas Esta sección esta obtenida de [MaMy].Suponemos una reacción química de forma que dos reactivos A y B dan lugar a unos productos C y D. dt dt A modo de ejemplo. en las cuales dos moléculas chocan para dar lugar a los productos.20338 kilogramos. podemos escribir que − d[NO] d[NO2 ] = = k[NO2 ][O2 ]. consideremos la reacción elemental A + B → C + D. 1. En general no puede establecerse una relación directa entre dichos números y los coeﬁcientes estequeométricos a y b.
Introducción a las ecuaciones diferenciales las concentraciones de los reactivos son 4 − x(t) y 2 − x(t).5. 7 Supongamos que tenemos dos recipientes conteniendo ambos una cierta sustancia en disolución. tenemos que c= Por otra parte. 14 . de donde y por lo tanto x(2) = 4 1− 1−2 1−2 ¡ 3 ¢2 2 1. tenemos que 1 1 log 3 = k + log 2. Consideremos el siguiente ejemplo. la ecuación d[C] = k[A][B]. Por otra parte. Los recipientes estan conectados entre sí. 2 Ã µ ¶! t 3 3 x(t) − 4 . permitiendo que entren sustancias del exterior y dejando salir también sustancias de los recipientes al exterior.1 Modelos de la química descritos por un sistema de ecuaciones Problemas de mezclas con varios recipientes 2 ¡ 3 ¢2 = 10 moles/litro. Podemos pensar por ejemplo en agua salada. como x(1) = 1. Además.5 1. 2 x(t) − 2 (4 − x(t))(2 − x(t)) y de la condición inicial x(0) = 0. de manera que puede pasar cierta cantidad de sustancia de un recipiente al otro y viceversa. 2 2 por lo que la constante de la reacción es k= Entonces 3 1 log . cada recipiente puede estar en contacto con el exterior. = t log + log 2 = log 2 log x(t) − 2 2 2 x(t) = 4 1− ¡ 3 ¢t 2 2 ¡ 3 ¢t . dt se reduce a x0 = k(4 − x)(2 − x). respectivamente. En este tipo de problemas se trata de determinar la concentración de sustancia disuelta en cada recipiente. 2 2 1 log 2. Resolvemos la ecuación diferencial Z Z x0 (t)dt 1 x(t) − 4 log = = kdt = kt + c. siendo x(0) = 0.
determinar la cantidad de sal en cada instante. De estos problemas. Para el caso del tanque A se tiene que y(t) ve = 6 l/m · 2 Kg/l + 1 l/m · Kg/l 100 y x(t) x(t) x(t) vs = 4 l/m · Kg/l + 3 l/m · Kg/l = 7 l/m · Kg/l 100 100 100 de donde obtenemos la ecuación diferencial x0 (t) = − 7 1 x(t) + y(t) + 12. Recordemos los problemas de mezclas con un único recipiente vistos anteriormente. cada uno con 100 litros de líquido se encuentran interconectados por medio de tubos. Para resolver el problema. La solucición (diluida) ﬂuye hacia el exterior del tanque A a razón de 4 l/m del tanque B a 2 l/m. 100 100 100 15 x(t) Kg/l 100 . Una solución de salmuera con una concentración de 2 Kg/l ﬂuye del exterior hacia el tanque A a razón de 6 l/m. El líquido contenido en el interior de cada tanque se mantiene bien agitado. Si inicialmente el tanque A contenía agua pura y el B 200 kg de sal.Introducción a las ecuaciones diferenciales Example 3 Dos grandes tanques. respectivamente. vemos que la variación de la cantidad de sal en A es x0 (t) = ve − vs donde ve es la velocidad de entrada de sal y vs es la velocidad de salida. 100 100 Procediendo de igual manera con el tanque B se tiene que y 0 (t) = ve − vs donde ahora ve = 3 l/m · y vs = 1 l/m · y(t) y(t) y(t) Kg/l + 2 l/m · Kg/l = 3 l/m · Kg/l. llamemos x(t) e y(t) las cantidades de sal en cada instante en los tanques A y B. El líquido ﬂuye del tanque A (ver dibujo posterior) hacia el tanque B a razón de 3 l/m y de B hacia A a razón de 1 l/m.
1. respectivamente. Otros problemas de este tipo se propondrán al ﬁnal del tema. 2 horas entre las zonas A y B y 5 horas entre la zona B y el exterior. 4 2 donde 1 U(t) = o C/1000 Kcal · 80. Example 4 Un ediﬁcio consta de dos zonas A y B (veáse la siguiente ﬁgura). ⎪ y (t) = ⎪ 100 100 ⎪ ⎩ x(0) = 0. 4 16 . = 3 3 y(t) 0 y 0 (t) − 100 100 Dejamos la resolución del sistema al alumno. y puede escribirse de forma matricial como µ 0 ¶µ ¶ µ ¶ ¶ µ 7 1 x (t) x(t) 12 − 100 100 + . y por consiguiente el problema de condiciones iniciales ⎧ ⎪ 0 ⎪ x (t) = − 7 x(t) + 1 y(t) + 12. Si la temperatura exterior es de 0 o C. como el del siguiente ejemplo. determinar la temperatura de cada zona. 100 100 El sistema es no autónomo.000 Kcal/h = 20 o C/h. llamemos x(t) e y(t) a las temperaturas de las zonas A y B. y(0) = 200. ⎪ ⎨ 100 100 3 3 0 x(t) − y(t).Introducción a las ecuaciones diferenciales de donde obtenemos la ecuación y 0 (t) = 3 3 x(t) − y(t). Esta aplicación se basaba en la ley de enfriamiento de los cuerpos de Newton. Vamos a ver qué pasa si el ediﬁcio tiene más de una estancia. Para resolver el problema.2 Climatización de ediﬁcios con varias estancias Anteriormente vimos una aplicación de las ecuaciones diferenciales a la climatización de ediﬁcios con una sola estancia.5. La zona A es calentada por un calefactor que genera 80000 Kcal/h. Las constantes de tiempo de transferencia de calor son entre la zona A y el esterior 4 horas. La capacidad caloríﬁca de la zona A es de 1/4 o C por cada 1000 Kcal. Entonces 1 1 x0 (t) = (0 − x(t)) + (y(t) − x(t)) + U(t).
1. H + + HOBr + Br− → H2 O + Br2 . Consideremos ahora la reacción hipotética 2A + B → 2C + D.3 Cinética de las reacciones químicas No todas las reacciones químicas son elementales. Por ejemplo. por lo que el proceso sigue las dos estapas siguientes: Br− + H + + H2 O2 → HOBr + H2 O. Otras reacciones pueden implicar muchos más pasos.Introducción a las ecuaciones diferenciales de donde conseguimos la ecuación 3 1 x0 (t) = − x(t) + y(t) + 20. dt dt dt dt donde los subíndices indican la variación de cada substancia como consecuencia de la primera reacción. la reacción en solución acuosa H2 O2 + 2Br− + 2H + → Br2 + 2H2 O no es el resultado del choque simultáneo de dos iones hidrógeno. A + E → D + C. en ninguna de las cuales se produce la interacción de más de tres partículas. 2 Dejamos la resolución del sistema al alumno. dos iones bromuro y una molécula de peróxido de hidrógeno (la probabilidad de que cinco especies estén en el mismo lugar al mismo tiempo es muy pequeña). de la que suponemos el siguiente mecanismo A + B → E + C. para la zona B tenemos 1 1 y 0 (t) = (0 − y(t)) + (x(t) − y(t)). 5 2 con lo que tenemos el sistema ½ x0 (t) = − 3 x(t) + 1 y(t) + 20. Por otra parte. dt dt dt dt 17 . de la segunda reacción − d[A2 ] d[E2 ] d[C2 ] d[D] =− = = = k2 [A][E]. 4 2 7 y 0 (t) = 1 x(t) − 10 y(t). Como sabemos de la primera reacción − d[B] d[C1 ] d[E1 ] d[A1 ] =− = = = k1 [A][B]. 4 2 Por otra parte. Otros ejemplos de este tipo serán estudiados en los ejercicios. aunque la mayoría de éstas pueden describirse mediante varias de ellas.5.
dt d[E] = k1 [A][B] − k2 [A][E] = [A](k1 [B] − k2 [E]). Este modelo da lugar a reacciones en las que no es predecible el comportamiento a largo plazo y que podríamos denominar “caóticos”.04 miligramos en un semana. dt dt dt d[E1 ] d[E2 ] d[E] = − . Y + D. 1. tenemos un sistema de dos ecuaciones diferenciales. dt d[Y ] = k2 [B][X] − k3 [X]2 [Y ]. 3X. dt dt dt d[C] d[C1 ] d[C2 ] = + . Se trataría de un sistema que no podríamos resolver analíticamente. La reacción del Brusselator La reacción del Brusselator sigue el siguiente esquema A B+X 2X + Y X → → → → X. [B]. ¿cuánto Torio tendremos al cabo de tres semanas? ¿Cuánto tiempo tiene que transcurrir para que la cantidad de Torio se reduzca a la mitad? 18 . El isótopo radioactivo del Torio 234 se desintegra a una rapidez proporcional a la cantidad existente en ese instante de tiempo. dt siendo ki las cosntantes de reacción de las cuatro reacciones que deﬁnen la reacción global. dt dt dt Entonces − d[A] = k1 [A][B] + k2 [A][E] = [A](k1 [B] + k2 [E]). y siento las concentraciones iniciales de X e Y nulas (ver [McPo]). donde las concentraciones [A].Introducción a las ecuaciones diferenciales Es obvio que d[A] d[A1 ] d[A2 ] = + . aunque tendría solución. E.6 Ejercicios 1. Si 100 miligramos de este material se reducen a 82. dt Si suponemos que la concentración del elemento B permanece constante. [D] y [E] son constantes y las concentraciones [X] e [Y ] siguen la ley d[X] = k1 [A] − k2 [B][X] + k3 [X]2 [Y ] − k4 [X].
El líquido ﬂuye del tanque A hacia el tanque B a razón de 4 litros por minuto y del tanque B al tanque A a 1 litro por minuto. Entonces se vierte en el tanque agua que contiene 0. Hacia el tanque A entra del exterior agua a razón de 3 litros por minuto 19 . ¿en qué momento la temperatura interior del ediﬁcio será de 16 o C? 8. La temperatura corporal de cadáver era de 350 C en dicho momento.055% de CO2 . y al minuto se enfrió a 880 C en un cuarto que está a 200 C. De observaciones experimentales se sabe que la temperatura superﬁcial de un objeto cambia con una rapidez proporcional a la diferencia de temperatura del objeto y su entorno. Dos tanques que contienen cada uno 50 litros de líquido se ecuentran interconectados por medio de dos tubos. En un día caluroso con una temperatura exterior de 40 o C.Introducción a las ecuaciones diferenciales 2. Un tanque contiene inicialmente 1000 litros de solución salina que contiene 10 Kg. de forma que por el otro extremo de la galería sale la misma cantidad de aire. Supongamos que decidís matar al profesor de ecuaciones diferenciales. Una vez perpetrado el hecho. El aprovechamiento es de medio grado por cada 1000 kilocalorías y la constante de tiempo del ediﬁcio es de 3 horas. A mediodía se apaga el calentador y la gente regresa a casa. 9. Determinar la sal que habrá en el tanque después de 20 minutos. de sal por litro se vierte en el tanque a la razón de 10 l/ min mientras que simultaneamente. la solución bien mezclada sale del tanque a razón de 15 l/ min . mientras que el aire del exterior tiene un 0. Este hecho es conocido como la ley de enfriamiento de Newton. de sal. La temperatura exterior permanece constante a 12 o C durante el resto de la tarde. ¿Qué concentración de CO2 habrá al cabo de 20 minutos en la galería? 7. En una mañana de sábado. y se deja que la mezcla salga del tanque con la misma rapidez. ¿cuánto tiempo debe de transcurrir para que se enfrie hasta los 650 C? 3. ¿A que hora se produjo el horripilante y brutal suceso? 4. determinar la temperatura del taller a lo largo del día. Un taller mecánico sin calefacción ni aire acondicionado tiene una constante de tiempo de 2 horas.05 kilogramos de sal por litro a una velocidad de 8 litros por minuto. 6. Encontrar la cantidad de sal que hay en el tanque en un momento t. Si inicialmente la temperatura del ediﬁcio era de 35 o C.2% de CO2 . Una hora más tarde la temperatura era de 330 C. Si la temperatura de una taza de café es de 950 C recién servida. Otra solución salina que contiene 25 Kg. ¿Cuál es el valor máximo de temperatura que puede tener el ediﬁcio en estas condiciones? 10. Un tanque contiene originalmente 400 litros de agua limpia. 5. mientras las personas trabajan. Si la temperatura exterior varía según la función E(t) = 30 − 15 cos(2πt/24). se enciende dentro de un ediﬁcio un aparato aire acondicionado que disipa 24000 kilocalorías por hora. determinar la temperatura al cabo de 3 horas. se encuentra el cuerpo en el despacho del mismo que está a una temperatura de 200 C a las 6 de la tarde. Se instalan ventiladores que introducen en la galería 9 metros cúbicos de aire del exterior por minuto. Si la constante de tiempo del ediﬁcio es de 3 horas. un calefactor mantiene la temperatura interior de un ediﬁcio a 21 o C.2m hay un 0. El líquido contenido en cada tanque se mantiene prefectamente agitado. En una galería subterranea de 15 × 15 × 1.
12. Una salmuera con concentración de 3 kilos por litro ﬂuye del exterior hacia el tanque A a razón de 5 litros por minuto. Dos grandes tanques. determinar la cantidad de sal en cada instante. 20 . Si inicialmente el tanque A contiene 25 kilos de sal y el tanque B no contiene nada de sal.Introducción a las ecuaciones diferenciales y la solución ﬂuye hacia el exterior por el tanque B a la misma velocidad. Solamente la zona A es calentada por un calefactor.000 kilocalorias por hora. saliendo hacia el exterior a la misma velocidad por un tubo situado en el tanque B. El líquido ﬂuye del tanque A hacia el B a razón de 5 litros por minuto. El líquido contenido en el interior de cada tanque se mantiene bien agitado. Las constantes de transferencia de calor son 4 horas entre la zona A y el exterior. determinar la cantidad de sal en cada instante de tiempo. Un ediﬁcio consta de dos zonas A y B. 11. cada uno de 50 litros se encuentran interconectados por un tubo. 5 horas entre la zona B y el exterior y 3 horas entre las dos zonas. ¿a qué temperatura puede llegar a enfriarse la zona B? Nota: las constantes de transferencia de calor son las inversas de las constantes que aparecen el la ley de enfriamiento de Newton. La capacidad caloríﬁca de la zona A es de 1/4 de grado Celsius por cada 1000 kilocalorias. Si la temperatura exterior es de 0 grados Centígrados. Si el tanque A contiene inicialmente 50 kilos de sal y el tanque B contiene 100 kilos. que genera 80.
¿a qué temperatura puede llegar a calentarse la zona del ático? 21 . El área habitacional es refrigerada por medio de una unidad de aire acondicionado de 2 toneladas que disipa 24000 kilocalorias por hora. 4 horas entre la zona B y el exterior y 4 horas entre ambas zonas. Para ﬁnes de refrigeración una casa consta de dos zonas: la zona de ático A y la zona B o habitacional. Si la temperatura exterior permanece a 40 grado centígrados. La constantes de transferencia de calor son 2 horas entre la zona A y el exterior.Introducción a las ecuaciones diferenciales 13. La capacidad caloríﬁca de la zona B es de 1/2 grado centígrado por cada 1000 kilocalorias.
Introducción a las ecuaciones diferenciales 22 .
pero son desconocidas. Basta considerar el problema ½ 0 2 y = ey ..1 Introducción (2. Consideramos un problema de condiciones iniciales de la forma ½ 0 y = f(t. t0 . y0 ). y(0) = 4.. y0 ) es suﬁcientemente diferenciable en un entorno de t0 . Sin embargo. Métodos de Taylor. i = 1. Los métodos que vamos a estudiar consisten en generar una sucesión y0 .1).. Es por ello importante determinar métodos que permitan obtener aproximaciones de dichas soluciones. Por ejemplo. y0 )h + y00 (t0 . y0 ) + y0 (t0 . Si t1 está en dicho entorno y denotando h = t1 − t0 . ..1) donde la función f : Ω ⊆ Rm+1 → Rm es suﬁciente regular para que dicho problema tenga solución ∂f única. Vamos a ver en este tema varias maneras de construir dicha sucesión. que existen. y). m continuas. para un cierto tf > t0 (análogamente se haría para tf < t0 ). y0 ) y buscamos cómo aproximar el valor de y(tf . y0 )h2 + . denotemos su solución por y(t. n! 23 . En esencia.. Análisis del error global. + 1! 2! 1 n) + y (t0 . . t0 .2 Métodos de Taylor Este método se basa en suponer que la solución y(t. Métodos de Runge—Kutta: tablas de Butcher. dada un problema de condiciones iniciales arbitrario.. t0 . y0 ) = y(t0 . Extrapolación de Richardson. f y ∂yi . y1 .Capítulo 2 Métodos de un paso Sumario. t0 . 2. t0 .. 2. entonces 1 1 y(t1 . es muy posible que no sepamos cómo hallar dicha solución. t0 . t0 . t0 .. yn de manera que yn sea un valor aproximado de y(tf . y0 )hn + O(hn ). dado el problema (2. y(t0 ) = y0 . y0 ). t0 .
y0 ) + fn−1 (t0 .. y0 ) = ∂t ∂y ∂ ∂ f2 (t0 . y0 ). y0 ) ∂y 1 0 denotamos el gradiente de f1 (t0 . 1! 2! n! Veamos qué forma particular tiene esta aproximación para diferentes valores de n. y0 ). y00 (t0 . y0 ) = fn−1 (t0 . esto es y1 = y0 + 1 1 1 f1 (t0 .. t0 . y0 ) = ∂t ∂y = f3 (t0 . y0 ) = f2 (t0 . y0 ) = y0 . y0 )hn 1! 2! n! es una aproximación de y(t1 .. Inductivamente. y0 )f(t0 . y0 ) + f2 (t0 . y0 ). y0 ). y0 )y0 (t0 .. y0 ) = d 0 d y (t0 ...2) . t0 . y0 ). t0 . + fn (t0 . y0 ) + fn−1 (t0 . y0 ) = ∂t ∂y = f2 (t0 . y(t1 . t0 . y0 )h2 + . y0 ) = d 00 d y (t0 . y0 ) = d n−1) d (t0 . y0 ) = fn−1 (t0 . entonces yn) (t0 . si yn−1) (t0 . y0 ) + f1 (t0 . y0 ) dt dt ∂ ∂ f2 (t0 . t0 . y0 ) = ∂t ∂y = fn (t0 . y0 ) + f1 (t0 .Métodos de un paso donde O(hn ) es denota una función g(h) para la cual existe una constante positiva k tal que |g(h)| ≤ k|hn |. t0 . y0 ) y dt dt ∂ ∂ fn−1 (t0 . t0 . y0 ) = ∂t ∂y ∂ ∂ fn−1 (t0 . y0 )hn + O(hn ). sustituyendo en la fórmula original y(t1 . y0 ) dt dt ∂ ∂ f1 (t0 . y0 )y0 (t0 . Así. donde por ∂ f (t . y0 )h2 + . y0 )h + f2 (t0 . y0 ) ≈ y1 = y0 + 24 (2. t0 . y0 )h + f2 (t0 . y0 )h + f2 (t0 . y0 ) = f1 (t0 . y0 (t0 . y0 )y0 (t0 . y0 )f(t0 . y0 ) = f(t0 . y0 )h2 + . y0 ) = ∂t ∂y ∂ ∂ f1 (t0 . + 1! 2! 1 fn (t0 . y0 )hn . Entonces y(t0 . t0 . t0 . + fn (t0 . t0 . y0 ). n! 1 1 1 f1 (t0 . t0 . y0 ). y0 ) = y0 + + con lo que 1 1 f1 (t0 . t0 . t0 . y3) (t0 . y0 ) + f2 (t0 . y0 )f(t0 . t0 . y0 ) = f1 (t0 .
es decir. y0 )(t − t0 ).2. que es la expresión (2. y0 ) ≈ y1 = y0 + f(t0 . t0 . t0 . Imaginemos que m = 1. t0 . Como vemos. Si sustituimos t por t1 en la recta anterior obtenemos y(t1 . y0 )h. 25 .Métodos de un paso Figura~2. t0 . t0 . 1! y tiene un claro signiﬁcado geométrico. 2. Veamos cómo funciona el método de Euler con un ejemplo.3) para ecuaciones de dimensión uno. y(0) = 1. y0 )(t − t0 ). y0 )|.1: El método de Euler.2) queda de la forma 1 (2. Consideremos el problema de condiciones iniciales ½ 0 y = y. y sustituyendo cada elemento de la expresión anterior por su valor obtenemos y − y0 = f (t0 . El error e1 es |y1 − y(t1 .1 nos muestra gráﬁcamente el método.1 Método de Euler El método de Taylor con n = 1. y0 ) ≈ y1 = y0 + f (t0 .3) y(t1 . t0 . y0 ) = y 0 (t0 . se trata de una ecuación diferencial. recibe el nombre de método de Euler y fue quizás el primer método numérico generado mucho antes de la existencia de ordenadores. la expresión (2. Entonces la recta tangente de la solución y(t. y0 ) para t = t0 tiene la forma y − y(t0 . y0 )h. La ﬁgura 2.
Para construir y2 . y7 + f (7h.1 · 0..1 = 1.1 = 1. y0 + f (0. 0.1 y calculamos y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 = = = = = = = = = = yi = yi−1 + f(ti−1 . 1) − y1 | = |e0. tiene por solución y(t.1 = 2.59374.1.61051 · 0. .. y6 )h = 1. n.14359 · 0. Tomamos una partición P del intervalo [t0 . entonces y1 = y0 + y0 h = 1 + 0. cantidad que suele llamarse tamaño de paso y n el número de pasos. el error cometido Si ahora.35795. Para reducir el error se procede de la siguiente manera. dado por y2 = y1 + f(t1 . 1). entonces h = 1 e y el error esto es.1 = 1. e1 = |y(1.61051. 1) − y1 | = |e − 2| ≈ 0. y9 + f (9h. y8 + f (8h.1 = 1. y2 )h = 1. < tn−1 < tn = tf . Construimos la sucesión yn de la siguiente manera y1 = y0 + f(t0 .1 + 1.35795 + 2.61051 + 1. y estimemos por el método de Euler y(0. 0. 26 ..77156 + 1. i = 0.77156. Esto se debe a que estamos tomando aproximaciones locales. 1. tomamos la aproximación mediante el método de Euler del problema ½ 0 y = f(t. suele tomarse hi = h.1. y8 )h = 2.00517092. y4 )h = 1..1.14359. y1 + f (h. yi−1 )h.. yi−1 )h = yi−1 + f(t0 + (i − 1)h.331 · 0. Ahora bien. 1) = et . En este caso el método de Euler queda como para i = 1. Tomemos t1 = 0. e1 = |y(0. y7 )h = 1.94872 · 0.1 − 1.1 = 1.1 = 1..21. y6 + f (6h.1 = 2. y5 )h = 1.718282. . En general..94872 + 1..331. y0 )h = 1 + 1 · 0. n − 1.4641 · 0. y1 es una aproximación de y(t1 . 0. y9 )h = 2.1. Como h = 0. .77156 · 0.. tomamos t1 = 1.1 = 1.4641 + 1.1 = 1. y3 + f (3h. n..331 + 1. En el ejemplo anterior. y3 )h = 1.1 = 2. el error aumenta considerablemente..1.35795 · 0. tomamos h = 0. y0 )h0 .. y de forma recurrente para i = 1.21 · 0. y1 )h = 1. . y5 + f (5h. y0 ). y1 = y0 + y0 h = 1 + 1 = 2. yi−1 )hi−1 . y(t1 ) = y1 .1. y1 )h1 . Como vemos.4641. Deﬁnimos hi = ti+1 − ti . i = 0.1| ≈ 0. t0 . 1. esto es P = t0 < t1 < t2 < . y2 + f (2h.21 + 1.14359 + 2.94872. n − 1. y4 + f (4h. 0.Métodos de un paso que como sabemos... tf ]. yi = yi−1 + f(ti−1 . y).
2.. y(0) = 1.051 0. y0 ). y0 ) f(t0 .2 Método de Taylor de orden 2 1 1 f1 (t0 . y0 ) = f(t0 . y0 ) ≈ y1 = y0 + f(t0 .718 0. t0 . al dividir el tamaño de paso h por diez.125 0. 2. yi−1 ) h2 . ∂t ∂y Si dividimos el intervalo [t0 .082 0. y0 ) + f(t0 .Métodos de un paso y ahora los errores son ei = |ei∗0... 1! 2! Si hacemos n = 2. de manera que estos van creciendo. 2 ∂t ∂y para i = 1. que puede estimarse con la recurrencia yi = yi−1 + f(t0 + (i − 1)h. Vemos no obstante que los errores se van acumulando desde e1 hasta e10 .102 0.01 h = 0. si disminuimos el tamaño de paso. para i = 1. Sin embargo.000014 Como vemos. el error ha disminuido notablemente. y0 ) + f(t0 .125 Como vemos. y0 )h + 2 µ ¶ ∂ ∂ f(t0 . y0 ) f(t0 . .013 0. 27 . yi−1 ) + f(t0 + (i − 1)h. Si consideramos el problema de condiciones iniciales anterior ½ 0 y = y. y0 )h2 .019 0.011 0.001 0.038 0.0001 h = 0. n. 10.005 0. yi−1 )h ¶ µ 1 ∂ ∂ + f(t0 + (i − 1)h.. que nos da la siguiente tabla aproximada e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 0. vemos en la siguiente tabla como los errores al ﬁnal dismuyen h = 1 h = 0. y0 ) h2 . tf ] en n intervalos igualmente espaciados siendo el tamaño de paso h = (tf − t0 )/n. a pesar de que en los pasos intermedios la aproximación del método de Euler no coincide en su condición inicial con la solución del problema de condiciones original.028 0. y0 ) y f2 (t0 . ∂t ∂y ∂ ∂ f(t0 . y0 ) con yn . y0 )h + f2 (t0 . y0 ) = se puede reescribir como 1 y(t1 . t0 . el valor y(tf .00014 0.. Veremos posteriormente una explicación a este hecho. yi−1 ) f(t0 + (i − 1)h. y0 ) ≈ y1 = y0 + y dado que f1 (t0 . el error ﬁnal aproximadamente también hace esta operación.1 − yi |.065 0..1 h = 0. t0 .00001 0. observamos que el método de Taylor queda de la siguiente forma: y(t1 .001 h = 0. .
0017 0.64745. para i = 1.5 · 10−9 4.001 h = 0.00017 0. y) = y.82043. y) ∂t =0y ∂f (t. al dividir por 10 el tamaño de paso.1 − yi |.. Además.0013 0.. por lo que de la forma = 1. 1+h+ y5 = 2 µ ¶ h2 y6 = 1+h+ y5 = 1.01 h = 0. Tomando h = 0.34923.0042 0. que nos da la siguiente tabla aproximada e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 0. se mejora notablemente el error con respecto al método de Taylor.5 · 10−11 Como vemos. 1+h+ y3 = 2 ¶ µ h2 y4 = y3 = 1.4909. 2 ¶ µ h2 y6 = 2.5 · 10−7 4. 2 ¶ µ h2 y8 = 2.45618. 28 .0009 0. 10. . los errores para diferentes tamaños de paso son h = 1 h = 0. siendo éste además de orden dos.01157.0006 0.0042 Como vemos.1 (n = 10) y calculando obtenemos µ ¶ h2 y1 = 1+h+ y0 = 1.Métodos de un paso tenemos que f (t. 2 ¶ µ h2 y1 = 1. y) ∂y para i = 1. el error es aproximadamente el del paso anterior al cuadrado. 1+h+ y9 = 2 ¶ µ h2 y9 = 2. n.2183 0.71408...22279.22103.1 h = 0.0028 0.00001 0. 1+h+ y2 = 2 ¶ µ h2 y2 = 1. el error decrece notablemente en comparación al obtenido al aplicar el método de Euler.105.0001 h = 0. 1+h+ y7 = 2 µ ¶ h2 y8 = 1+h+ y7 = 2. 2 2 ∂f (t. .. 1+h+ 2 ¶ µ h2 y4 = 1.0022 0.0004 0.0034 0.. es decir. y así la recurrencia anterior se expresa ¶ µ 1 h2 2 yi = yi−1 + yi−1 h + yi−1 h = 1 + h + yi−1 . 1+h+ y10 = 2 y ahora los errores son ei = |ei∗0.000045 4.
y(t0 )) + f(tf . t0 . y0 ) = y(t. t0 . seguimos disminuyendo el error producido. t0 . t0 . y(s))ds ≈ [f(t0 . Entonces y1 = y0 + h ∗ [f(t0 . 2 t0 siendo h = tf −t0 . y0 ) es solución y(t. y(s))ds t0 mediante algún método de integración numérica apropiado. t0 . y0 ) es conocido. y que se conocen como métodos de Runge—Kutta. obtenemos éste por un algoritmo de los estudiados anteriormente para tamaño de paso h. y0 ). y0 ) es precisamente el valor que tenemos que aproximar mediante el método. y0 )ds = f(s. y). y1 )] . y(tf ))] . y0 ) + f(tf . t0 . y0 ) − y0 Z t Z t = y(s. y esto frecuentemente no es tarea fácil. Sin embargo f(tf . y(s))ds. A modo de primer ejemplo. y(tf )) es desconocido dado que y(tf ) = y(tf . supongamos que dicha integral se aproxima mediante el método del trapecio. Además. que se construye de la siguiente manera. y0 ))ds. el método de Taylor presenta el problema de que hay que derivar sucesivamente las funciones que determinan la ecuación o sistema de ecuaciones diferenciales. Como y(t. que es el método de Euler. Sin embargo. y(t0 )) = f(t0 . 2. t0 t0 Entonces. y(s. el incremento del orden no mejora el algoritmo dado que fm (t. esto es Z tf h f(s. t0 . Veamos a continuación una familia de métodos que presentan un avance en este sentido. en el caso del ejemplo anterior para la ecuación y 0 = y. El valor f(t0 . Para obtener un valor aproximado de dicho valor para poder aplicar el método. y(t) − y0 = t0 Los métodos de Runge—Kutta se basan en obtener aproximaciones de la integral Z t f(s. y) = 0 para m ≥ 3. y(t0 ) = y0 . 2 29 .Métodos de un paso Aumentando el orden del método de Taylor. por ejemplo ∗ y1 = y0 + hf(t0 . dicha función puede calcularse a partir de la ecuación integral Z t f(s.3 Métodos de Runge—Kutta de orden dos Este método se basa en la ecuación integral asociada al problema de condiciones iniciales ½ 0 y = f(t. y0 ) − y(t0 .
y0 ) que buscábamos.0009 0. 2 ½ ∗ y9 = (1 + h) y8 = 2.34923. ∗ y8 = y7 + h (y7 + y8 ) = 2. ∗ y4 = y3 + h (y3 + y4 ) = 1.48416.0028 0. y(0) = 1. yi ) + f(t0 + (i + 1)h. se dice que es un método de Runge—Kutta de dos etapas. ∗ y6 = y5 + h (y5 + y6 ) = 1.0013 0.0042 30 .45618. ⎩ yi = yi−1 + 1 (g1 + g2 ). ∗ y7 = y6 + h (y6 + y7 ) = 2.34313.4909. se tiene que de forma compacta ∗ yi+1 = y0 + hf(t0 + ih. 2 ½ ∗ y6 = (1 + h) y5 = 1.00247. yi+1 ) yi+1 = y0 + 2 que se conoce como método de Heun. Por ello. ∗ y5 = y4 + h (y4 + y5 ) = 1. hay dos etapas.0006 0. 2 cuyos errores son e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 0.Métodos de un paso será la aproximación de y(tf .21273. 2 ½ ∗ y2 = (1 + h) y1 = 1. 2 ½ ∗ y3 = (1 + h) y2 = 1. yi−1 + g1 ).105.82043.0004 0. ∗ y9 = y8 + h (y8 + y9 ) = 2. 2 ½ ∗ y4 = (1 + h) y3 = 1. ¤ h£ ∗ f(t0 + ih.1. yi−1 ). 2 Veamos cómo se implementa este método en nuestro ejemplo de costumbre ½ 0 y = y.44507. Suele escribirse de forma más compacta como ⎧ ⎨ g1 = hf(ti−1 .7018. 2 ½ ∗ y7 = (1 + h) y6 = 2. yi ).0022 0.81219.71408. una inicial donde se calcula ∗ yi y otra posterior donde ya se obtiene la aproximación propiamente dicha. ∗ y10 = y9 + h (y9 + y10 ) = 2. Si tomamos un tamaño de paso h = (tf −t0 )/n. ∗ y1 = y0 + h (y0 + y1 ) = 1.2155. 2 ½ ∗ y8 = (1 + h) y7 = 2. 2 ½ ∗ y10 = (1 + h) y9 = 2.0034 0. Como vemos.22279.0017 0.22103.64745. ∗ y3 = y2 + h (y2 + y3 ) = 1. y como veremos posteriormente de orden dos. 2 ½ ∗ y5 = (1 + h) y4 = 1. t0 .00017 0.1 son ½ ∗ y1 = (1 + h) y0 = 1. g2 = hf(ti−1 + h.63999. ∗ y2 = y1 + h (y1 + y2 ) = 1. Los valores que obtenemos para tamaño de paso h = 0.01157.
2 2 e h 1 yi = yi−1 + hg2 = yi−1 + hf(ti−1 + . 2 ∂t ∂y e igualando coeﬁcientes obtenemos que ⎧ ⎨ b1 + b2 = 1. yi−1 )h + f(ti−1 . La forma más general posible para un método de Runge—Kutta de orden dos es ⎧ ⎨ g1 = hf(ti−1 .Métodos de un paso Obsérvese que son similares a los obtenidos en el método de Taylor de segundo orden. yi−1 + 1 g1 ). yi−1 + g1 ) 2 2 h h = yi−1 + hf(ti−1 + .1 h = 0. Si variamos el tamaño de paso.000045 4. yi−1 + f(ti−1 . yi−1 ) + a21 (ti−1 . ti−1 . ⎩ b2 a21 = 1/2. yi−1 ) .00001 0. yi−1 ) h2 .5 · 10−9 4. yi−1 ) + f(ti−1 . la aproximación mediante la serie de Taylor de orden dos de y(t.01 h = 0. yi−1 ). yi−1 ) + b2 h c2 (ti−1 . yi−1 + a21 g1 ). Así. ⎩ yi = yi−1 + b1 g1 + b2 g2 . lo cual nos proporciona una familia de métodos de Runge—Kutta de orden dos según los valores de b2 . y b2 6= 0. b2 c2 = 1/2.5 · 10−7 4. tenemos ½ g1 = hf(ti−1 . g2 = hf(ti−1 + h . obtenemos el método de Heun anteriormente descrito. yi−1 )f(ti−1 . 2 2 31 .001 h = 0. yi−1 ) + O(h). ∂t ∂y por lo que yi = yi−1 + b1 g1 + b2 g2 µ ¶ ∂f ∂f = yi−1 + (b1 + b2 )hf(ti−1 .2183 0.0001 h = 0. g2 = hf(ti−1 + c2 h. donde los tiempos ti−1 no tienen porqué ser uniformemente distribuidos. cuando b2 = 1/2. yi−1 ) era ¶ µ 1 ∂ ∂ yi = yi−1 + f(ti−1 . 1 Como b2 6= 0. ∂t ∂y 2 Por otra parte. yi−1 ). obtenemos los siguientes errores para los siguientes valores h = 1 h = 0. Cuando b2 = 1.5 · 10−11 que reproducen los obtenidos en el método de Taylor anteriormente mencionado.0042 0. Si tomamos g2 (h)/h y desarrollamos mediante la serie de Taylor de primer orden obtenemos 0 g2 (h) = g2 (0) + hg2 (0) + O(h) µ ¶ ∂f ∂f = f(ti−1 . yi−1 )). yi−1 )f(ti−1 . tenemos que a21 = c2 = 2b2 y b1 = 1 − b2 . yi−1 ) + h c2 (ti−1 . yi−1 ) f(ti−1 . yi−1 ) + a21 (ti−1 .
j = 2. + amm−1 gm−1 ). . b = (b1 .... y(0) = 1. y como veremos posteriormente... . bm ) y A = (ajk ) ∈ Mm×m (R) con ajk = 0 si k > j. j = 1...... m y ajk ... yi−1 + 1 g1 ).. que es el método clásico de Runge—Kutta de cuatro etapas. ⎨ g3 = hf(ti−1 + c3 h.. = hf(ti−1 + h ... un método de Runge—Kutta explícito de m etapas es de la forma yi = yi−1 + donde m X j=1 bj gj . j = 1. a32 . yi−1 ). yi−1 )... m. . . 2 2 = hf(ti−1 + h . cm y en forma matricial ct A b a21 a31 . bj . ⎪ ⎪ ⎪ g2 = hf(ti−1 + c2 h. ⎪ ⎪ . 32 que da ....... am1 b1 ⎧ ⎪ g1 = hf(ti−1 .. yi−1 + am1 g1 + am2 g2 + . .. estos coeﬁcientes se agrupan según la tabla 0 c2 c3 . los coeﬁcientes del método. Si aplicamos este método a nuestro ejemplo ½ 0 y = y. . Así. siendo cj . bm−1 bm donde c = (0...... yi−1 + 1 g2 ). el método dado por la tabla 0 1 2 1 2 1 2 1 y concretado en ⎧ ⎪ g1 ⎪ ⎨ g2 ⎪ g3 ⎪ ⎩ g4 0 0 1 6 1 2 0 1 3 1 1 3 1 6 = hf(ti−1 . yi−1 + a21 g1 ). En general.. .. cuarto orden.. c2 . 2 2 = hf(ti−1 + h...... yi−1 + g3 ).. m.. am2 b2 . b2 ... ⎪ ⎪ ⎩ gm = hf(ti−1 + cm h.. 1 1 1 1 yi = yi−1 + g1 + g2 + g3 + g4 6 3 3 6 es el método de Kutta de 1905.. j − 1.Métodos de un paso que es el algoritmo de Runge de 1895. amm−1 . Normalmente. k = 1. cm )... yi−1 + a31 g1 + a32 g2 ).
t0 .82212.2214.0001 6. y1 .00001 6.1 · 10−6 h = 0.4 Análisis del error en los métodos de orden uno Consideramos un problema de condiciones iniciales de la forma ½ 0 y = f(t. y).. = 2. En general.1 son y1 y3 y5 y7 y9 cuyos errores son e1 e3 e5 e7 e9 = 8. Este hecho sí se aprecia en los tamaño de paso hasta 0. ﬁjado t1 > t0 y un tamaño de paso h = t1 −t0 . y8 = 2... los métodos que conocemos son de la forma yi = yi−1 + hΦ(ti−1 .1 · 10−7 = 6..001 1. • Errores de redondeo al trabajar los computadores con precesión ﬁnita... yn anteriormente mencionados. y10 = 2. donde la función f : Ω ⊆ Rm+1 → Rm es suﬁciente regular para que dicho problema tenga solución única. h). Como hemos visto hasta ahora. y(t0 ) = y0 .Métodos de un paso tenemos que los valores que obtenemos para tamaño de paso h = 0. construir una sucesión y0 .34986.4 · 10−8 = 3.00995 2. 33 .3 · 10−7 = 1.1 · 10−6 = 1. los errores de redondeo hacen que no se aprecie que el error del paso h/10 es aproximadamente el del paso h elevado a la cuarta potencia. dentro de los errores matemáticos podemos distinguir los siguientes tipos. yn de manera que n sean una aproximación de la solución y(t. obtenemos los siguientes errores para los siguientes valores h=1 h = 0. y6 = 1.22554.71828. = 1. = 2.64872. y4 = 1.. 2.01 2. n. = 1.3 · 10−7 e8 = 1.22 · 10−15 h = 0. y1 .7 · 10−6 e2 = 1.4596. Obsérvese que son similares a los obtenidos en el método de Taylor de segundo orden. yi−1 .26 · 10−14 Como vemos..49182. y2 = 1.4 · 10−6 e10 = 2.01375. Si variamos el tamaño de paso. y0 ) en los tiempos ti = t0 + hi.. i = 1. . los métodos numéricos de Taylor y Runge—Kutta se basan en. Como hemos puesto de maniﬁesto con algunos ejemplos. En general. (2.001.8 · 10−7 e4 = 4.25 · 10−10 h = 0. .4) que dan lugar a los valores y0 .10517. .6 · 10−7 e6 = 8.1 0. estos métodos tienen inherentemente asociados unos errores que se deben a dos causas bien diferenciadas: • Errores matemáticos debidos al método numérico empleado.38 · 10−14 h = 0.1 · 10−6 = 1.
Básicamente. 34 . yi−1 ). h) lim h→0 = −y0 (ti−1 . y0 ) + f(ti−1 . t0 . tienen gran importancia dos errores locales que a continuación describimos. dicho método es consistente si n X j=1 bj = 1.5) es Φ(t.. h) − y(ti . y el método numérico en cuestión se dirá convergente si h→0 0≤i≤n lim max ||ei || = 0. En el estudio del error global de un método. y0 ).4) se dice consistente si Φ(t. Se llama error local del paso i como li = yi − y(ti . El último tipo de error local que vamos a introducir es lo que llamaremos el error local de truncamiento. y0 ) − y(ti−1 + h. y0 ) = lim h→0 h h→0 h + lim Φ(ti−1 . t0 . y0 )) = 0. y0 ). es decir. ti−1 . t0 . Otro concepto importante es el de consistencia. t0 . i = 0. t0 . t0 . yi−1 . t0 . t0 . aquel error que se obtendría al sustituir la solución real del problema de condiciones iniciales y(t. y0 ).5) Por ejemplo. En el caso de los métodos de Runge—Kutta. n. y0 ). y). En efecto ti y(ti−1 . y(ti−1 ) = yi . y(ti−1 . es decir. y. la función dada por (2. y). 1. t0 .Métodos de un paso Deﬁnimos el error global de la solución aproximada como ei = yi − y(ti . (2. la convergencia implica que el error global tiende a cero cuando lo hace el tamaño de paso. h). y). y. Un método numérico dado por (2. que deﬁnimos como ti = y(ti−1 . y(ti−1 . la diferencia entre el valor proporcionado por el método yi y el valor exacto proporcionado por la solución exacta del problema de condiciones iniciales ½ 0 y = f(t. 0) = n X hj−1 j=1 j! yj) = y0 = f(t. t0 . y0 ) en el método numérico implementado yi = yi−1 + hΦ(ti−1 . Si el método numérico es consistente. en el caso de los métodos de Taylor de orden n. 0) = f(t. entonces el error de local de truncaminento converge a cero cuando se divide por h... . y0 ) + hΦ(ti−1 . y(ti−1 .
Métodos de un paso por ser y(ti−1 ; t0 , y0 ) solución de la ecuación diferencial. Básicamente, la consistencia indica que ||ti || es al menos O(h2 ). Diremos que un método es consistente de orden n si es consistente y el error local de truncamiento es de orden O(hn+1 ). Vamos a analizar este último tipo de error para los método numéricos que conocemos.
Error local de truncamiento en el método de Taylor
1 1 1 f1 (ti−1 , yi−1 )h + f2 (ti−1 , yi−1 )h2 + ... + fn (ti−1 , yi−1 )hn 1! 2! n! n X hj fj (ti−1 , yi−1 ). = yi−1 + j! j=1
Recordemos que el método de Taylor es de la forma yi = yi−1 +
Entonces, el error local de truncamiento es ti = y(ti−1 ; t0 , y0 ) + hΦ(ti−1 , y(ti−1 ; t0 , y0 ), h) − y(ti ; t0 , y0 ) n ∞ X hj X hj j) y (ti−1 ; t0 , y0 ) − yj) (ti−1 ; t0 , y0 ) = j! j! j=1 j=1
∞ X hj yj) (ti−1 ; t0 , y0 ) = − j! j=n+1
hn+1 n+1) (ti−1 + ηh; t0 , y0 ), y (n + 1)!
con η ∈ (0, 1). Entonces, existe A > 0 de manera que ||ti || ≤ Ahn+1 , o equivalentemente ||ti || = O(hn+1 ), con lo que el método de Taylor truncado en el paso n es de orden n + 1.
Error local de truncamiento en los métodos de Runge—Kutta
Por simplicidad, vamos a considerar sólamente los métodos de dos etapas dados por ⎧ ⎨ g1 = hf(ti−1 , yi−1 ), g2 = hf(ti−1 + c2 h, yi−1 + a21 g1 ), ⎩ yi = yi−1 + b1 g1 + b2 g2 , donde a21 = c2 =
1 2b2
y b1 = 1 − b2 . Entonces
Φ(t, y, h) = b1 f(t, y) + b2 f(t + c2 h, y + a21 hf(t, y)). 35
Métodos de un paso Ahora bien ti = y(ti−1 ; t0 , y0 ) + hΦ(ti−1 , y(ti−1 ; t0 , y0 ), h) − y(ti ; t0 , y0 ) ∞ X hj yj) (ti−1 ; t0 , y0 ) = y(ti−1 ; t0 , y0 ) + hΦ(ti−1 , y(ti−1 ; t0 , y0 ), h) − j! j=0 ! Ã ∞ j−1 Xh = h Φ(ti−1 , y(ti−1 ; t0 , y0 ), h) − yj) (ti−1 ; t0 , y0 ) j! j=1 = h3 3) y (ti−1 + ηh; t0 , y0 ), 3!
En general, el error local de truncamiento en un método de Runge—Kutta de n pasos es O(hn+1 ).
debido a los valores de los coeﬁcientes del método de Runge—Kutta, siendo η ∈ (0, 1). De esta manera, vemos que ||ti || = O(h3 ).
Relación entre el error local de truncamiento y el error global
Vamos a ver cómo podemos controlar el error global a partir del error local de truncamiento. Como sabemos ei = yi − y(ti ; t0 , y0 ). Entonces ei = yi−1 + hΦ(ti−1 , yi−1 , h) + ti − y(ti−1 ; t0 , y0 ) − hΦ(ti−1 , y(ti−1 ; t0 , y0 ), h) = ti + ei−1 + h (Φ(ti−1 , yi−1 , h) − Φ(ti−1 , y(ti−1 ; t0 , y0 ), h)) . Sea M > 0 tal que ||ti || < M y teniendo en cuenta que h > 0, obtenemos la acotación ||ei || ≤ M + ||ei−1 || + h||Φ(ti−1 , yi−1 , h) − Φ(ti−1 , y(ti−1 ; t0 , y0 ), h)||. (2.6) Supongamos ahora que Φ satisface una condición de Lipschitz con constante L > 0 en la variable y, esto es ||Φ(ti−1 , yi−1 , h) − Φ(ti−1 , y(ti−1 ; t0 , y0 ), h)|| ≤ L||yi−1 − y(ti−1 ; t0 , y0 )||, ||ei || ≤ M + ||ei−1 || + h||Φ(ti−1 , yi−1 , h) − Φ(ti−1 , y(ti−1 ; t0 , y0 ), h)|| ≤ M + ||ei−1 || + hL||yi−1 − y(ti−1 ; t0 , y0 )|| = M + ||ei−1 ||(1 + hL). Aplicando la anterior desigualdad para i ≥ 1 tenemos ||e1 || ≤ ||e0 ||(1 + hL) + M, ||e2 || ≤ ||e1 ||(1 + hL) + M ≤ ||e0 ||(1 + hL)2 + M(1 + (1 + hL)), 36
con lo que la expresión (2.6) se reduce a
Métodos de un paso y en general
i−1 X (1 + hL)j . ||ei || ≤ ||e0 ||(1 + hL) + M i j=0
i−1 X (1 + hL)i − 1 , (1 + hL)j = hL j=0
sustituyendo en la expresión anterior, (1 + hL)i − 1 ||ei || ≤ ||e0 ||(1 + hL)i + M µ ¶ hL M M − . = (1 + hL)i ||e0 || + hL hL De la noción de exponencial real tenemos que 0 ≤ (1 + hL)i ≤ ehLi , y teniendo en cuenta que e0 = 0, concluimos que ||ei || ≤ M hLi (e − 1). hL
Ahora bien ih = tf −t0 , donde tf es el tiempo ﬁnal donde deseamos conocer la solución de la ecuación diferencial, por lo que M L(tf −t0 ) (e − 1). ||ei || ≤ hL Como por otra parte, M = Ahn+1 donde A > 0 [M ≈ O(hn+1 )], tenemos que ||ei || ≤ Bhp , por lo que acabamos de probar el siguiente resultado: Theorem 2 Si ||ti || ≈ O(hn+1 ) entonces ||ei || ≈ O(hn ).
Relación entre el error local y el error local de truncamiento
Ambos errores son muy parecidos en magnitud. Como se puede comprobar li = yi − y(ti ; ti−1 , yi−1 ) = yi−1 + hΦ(ti−1 , yi−1 , h) − = h Φ(ti−1 , yi−1 , h) − Ã X hj
yj) (ti−1 ; ti−1 , yi−1 ) !
X hj−1
y (ti−1 ; ti−1 , yi−1 ) .
(2. tendremos µ ¶p h ei (h/2) = y(ti . si multiplicamos 2p ei (h/2) = 2p y(ti . ti−1 . h/2) − y(ti . h) + y(ti . h) es la solución aproximada del problema de condiciones iniciales ½ 0 y = f(t. h/2) − 2p y(ti .7) donde la función d(t) no depende de h. que el error global tanto en los métodos de Taylor como en los de Runge—Kutta tiene esta estructura. y0 ) = d(ti )hp + O(hp+1 ). y0 ) ! 2. y). Si calculamos ahora la aproximación con paso h/2. 2 con lo que µ ¶p h y(ti . y0 ). y(t0 ) = y0 . t0 . y(ti−1 . t = ti = t0 + ih. t0 . t0 .Métodos de un paso Entonces li − ti = h Φ(ti−1 .5 Extrapolación de Richardson Supongamos que y(ti ) = y(ti . h) − 2p y(ti . yi−1 . y0 ) = O(hp+1 ). h/2) − y(ti . esto es. p 1−2 2 Por otra parte. y0 ) − = O(hp+1 ). con tamaño de paso h. h/2) y(ti . t0 . t0 . (2. aunque queda lejos de los objetivos de este curso.8) 2 Restando ei (h/2) − ei (h) = y(ti . p 1−2 38 . h/2) − y(ti . t0 . h) − y(ti . Puede comprobarse. h) = d(ti ) + O(hp+1 ). h) − Ã Ã X hj−1 j≥1 j! yj) (ti−1 . y0 ) = d(ti ) + O(hp+1 ). h) − j! yj) (ti−1 . h/2) − 2p y(ti . y0 ) = d(ti )hp + O(hp+1 ). y0 ) − y(ti . de donde y(ti . Supongamos que el error global ei (h) = y(ti . t0 . yi−1 ) X hj−1 j≥1 ! −h Φ(ti−1 . y calculando 2p ei (h/2) − ei (h) = 2p y(ti . t0 . h) µ ¶ 1 p = d(ti )h − 1 + O(hp+1 ) p 2 µ ¶p h = d(ti ) (1 − 2p ) + O(hp+1 ).
. m. j = 2. h/2) = 2 · 2. m y ajk . t0 .. Esto es lo que se conoce como el método de extrapolación de Richardson. yi−1 ). .. j − 1. un método de Runge—Kutta explícito de m etapas es de la forma yi = yi−1 + donde m X j=1 bj gj ... bm−1 39 bm .. 1−2 que nos da un error de 6. 0. ⎪ ⎪ .. h) − 2y(ti .. am1 b1 ⎧ ⎪ g1 = hf(ti−1 .. am2 b2 . + amm−1 gm−1 ).. 1) por el método de Euler con tamaño de paso h = 10000 y h/2 = 5000. aumentamos en uno el orden de convergencia. y0 ) que tiene al menos orden O(hp+1 ). m. .. ⎪ ⎪ ⎩ gm = hf(ti−1 + cm h.71821 − 2. y calculemos y(1. y(t0 ) = y0 . j = 1... j = 1..22853 · 10−9 . .. Teniendo en cuenta que el método de Euler es de orden uno. De este modo.000135902 y e(h/2) = 0.. y(0) = 1..0000679539. a32 .. yi−1 + am1 g1 + am2 g2 + .. siendo cj . En general... estos coeﬁcientes se agrupan según la tabla 0 c2 c3 .71815 = 2.... calculamos y(ti . cm a21 a31 .71828. yi−1 + a21 g1 )... A modo de ejemplo. con lo que el error ha disminuido notablemente con una serie de operaciones sencillas... 2. ⎨ g3 = hf(ti−1 + c3 h. .. yi−1 + a31 g1 + a32 g2 )... ⎪ ⎪ ⎪ g2 = hf(ti−1 + c2 h. los coeﬁcientes del método. amm−1 .. h/2) 1 − 2p es una aproximación de y(ti . k = 1. tomemos el problema ½ 0 y = y... .. Obtenemos los errores e(h) = 0.. y).6 Más sobre los métodos Runge—Kutta Partimos de la ecuación diferencial con condiciones iniciales ½ 0 y = f(t. bj .. h) − 2p y(ti .. sin por ellos aumentar la cantidad de operaciones de una forma drástica... Normalmente..Métodos de un paso por lo que y(ti .
yi−1 ). ti−1 . Como sabemos. yi−1 ) + f(ti−1 . la aproximación mediante la serie de Taylor de orden tres de y(t. yi−1 )2 2 f(ti−1 . Se satisfacen en general las condiciones de simpliﬁcación cj = j−1 X k=1 cjk . yi−1 + a21 g1 ) + f(ti−1 + c2 h.Métodos de un paso y en forma matricial ct A b donde c = (0. b2 . yi−1 ) 6 ∂t2 ∂y∂t ∂t ∂y ¶ 2 ∂ ∂ (2. su tabla de Butcher será 0 c2 c3 y el método será de la forma g1 = hf(ti−1 . bm ) y A = (ajk ) ∈ Mm×m (R) con ajk = 0 si k > j. ∂y ∂y Tomamos la función G2 (h) = f(ti−1 + c2 h. yi−1 ) + 2f(ti−1 .. . g2 = hf(ti−1 + c2 h. cm ). yi−1 ) h3 + O(h4 ).11) + f(ti−1 .10) 2 ∂t ∂y µ 1 ∂2 ∂2 ∂ ∂ + f(ti−1 . yi−1 + a21 g1 )a21 f(ti−1 . (2. b = (b1 . yi−1 ) f(ti−1 . . yi−1 ) era µ ¶ 1 ∂ ∂ f(ti−1 ..1 El método de 3 etapas Veamos cómo se genera el método de Runge—Kutta de tres etapas que tendrá error local de truncamiento ti ≈ O(h4 ). yi−1 ) + f(ti−1 . m. yi−1 ) f(ti−1 . yi−1 ) h2 yi = yi−1 + f(ti−1 . yi−1 + a21 g1 )a21 f(ti−1 . yi−1 )h + (2. yi−1 + a21 g1 ). yi−1 ). yi−1 + a21 g1 ).. yi−1 )2 . 21 ∂y 40 .. 2. yi−1 + a21 g1 )a2 f(ti−1 .9) Por otra parte... yi−1 ) f(ti−1 .. g3 = hf(ti−1 + c3 h. c2 .6. yi−1 )2 + f(ti−1 . yi−1 + a21 g1 ) + 2c2 f(ti−1 + c2 h. y yi = yi−1 + b1 g1 + b2 g2 + b3 g3 ... j = 2. . yi−1 + a31 g1 + a32 g2 ). yi−1 ) f(ti−1 . ∂t ∂y a21 a31 b1 a32 b2 b3 ∂2 ∂2 f(ti−1 + c2 h. derivamos dos veces respecto de h G02 (h) = c2 G00 (h) = c2 2 2 ∂ ∂ f(ti−1 + c2 h. yi−1 ) ∂t2 ∂t∂y ∂2 + 2 f(ti−1 + c2 h.
yi−1 ) + a32 (G2 (h) + hG02 (h))) ∂y ∂ + f(ti−1 + c3 h. yi−1 ) 2 ∂y y por otra parte G3 (h) = f(ti−1 + c3 h. yi−1 ) + 2c3 (a31 + a32 ) ∂t2 ∂t∂y ∂2 + (a31 + a32 )2 2 f(ti−1 . yi−1 )f(ti−1 . ∂y ∂t ∂y por lo que el desarrollo de Taylor en cero es 1 G3 (h) = G3 (0) + G03 (0)h + G00 (0)h2 + O(h3 ) 3 µ 2 ¶ ∂ ∂ f(ti−1 . yi−1 ) h = f(ti−1 . yi−1 ) f(ti−1 . yi−1 ) + c2 f(ti−1 . yi−1 ) f(ti−1 . yi−1 ) + c3 f(ti−1 . yi−1 )f(ti−1 . yi−1 ) . yi−1 )2 ∂y µ ¶ ∂ ∂ ∂ +2a32 f(ti−1 . yi−1 + a31 g1 + a32 hG2 (h)) (a32 (2G02 (h) + hG00 (h))) . yi−1 ) h2 + O(h3 ) + a21 f(ti−1 . yi−1 ) f(ti−1 . yi−1 ) h ∂t ∂y µ 2 2 ∂ 1 2∂ c2 2 f(ti−1 . yi−1 ) + a21 f(ti−1 . yi−1 ) + 2c2 a21 f(ti−1 . yi−1 + a31 g1 + a32 hG2 (h)) · (a31 f(ti−1 . yi−1 ) + 2 ∂t ∂t∂y ¶ 2 2 2 ∂ f(ti−1 . yi−1 + a31 g1 + a32 hG2 (h)). yi−1 ). yi−1 + a31 g1 + a32 hG2 (h)) + ∂t2 ∂2 f(ti−1 + c3 h. yi−1 ) + a32 (G2 (h) + hG02 (h))) . 2 ∂y G03 (0) = c3 ∂ ∂ f(ti−1 . yi−1 ) c2 f(ti−1 . yi−1 ) + a32 (G2 (h) + hG02 (h))) +2c3 ∂t∂y ∂2 2 + 2 f(ti−1 + c3 h. yi−1 ) + (a31 + a32 ) ∂t ∂y 41 . yi−1 + a31 g1 + a32 hG2 (h)) · (a31 f(ti−1 . yi−1 + a31 g1 + a32 hG2 (h)) ∂t ∂ + f(ti−1 + c3 h. ∂y G00 (h) = c2 3 3 ∂2 f(ti−1 + c3 h. ∂t ∂y de donde G00 (0) = c2 3 3 ∂2 ∂2 f(ti−1 . yi−1 )f(ti−1 . yi−1 ) f(ti−1 .Métodos de un paso de donde el desarrollo de Taylor de orden dos es 1 G2 (h) = G2 (0) + G02 (0)h + G00 (0)h2 + O(h2 ) 2 µ 2 ¶ ∂ ∂ = f(ti−1 . yi−1 )f(ti−1 . yi−1 ) + a21 f(ti−1 . yi−1 ) + (a31 + a32 ) f(ti−1 . yi−1 + a31 g1 + a32 hG2 (h)) · (a31 f(ti−1 . que derivando dos veces nos da G03 (h) = c3 ∂ f(ti−1 + c3 h.
2 b2 c2 + b3 c2 = 1 . tenemos de las ecuaciones ½ b2 c2 + b3 c3 = 1 . 2 b2 a21 + b3 (a31 + a32 ) = 1 . 3 b2 a2 + b3 (a31 + a32 )2 = 1 . yi−1 ) c2 f(ti−1 . ⎪ ⎪ ⎨ c3 = a31 + a32 .Métodos de un paso µ 1 2 ∂2 ∂2 + c3 2 f(ti−1 . Si tomamos c2 y c3 como parámetros. 6 De las dos últimas obtenemos que c2 = a21 con lo que usando la segunda y tercera. 2 ⎪ ⎪ ⎪ b2 c2 + b3 c2 = 1 . 2 3 3 ¯ 1 ¯ 2 ¯ 1 ¯ b2 = ¯ 3 ¯ c2 ¯ 2 ¯ c2 c3 c2 3 c3 c2 3 ¯ ¯ ¯ ¡ ¢ ¯ c3 c23 − 1 3 ¯= ¯ c2 c3 (c3 − c2 ) ¯ ¯ 42 que . yi−1 )f(ti−1 . yi−1 ) 2 ∂t ∂t∂y ∂2 + (a31 + a32 )2 2 f(ti−1 . 6 b3 a32 a21 = 1 . yi−1 )2 ∂y µ µ ¶¶ ∂ ∂ ∂ + 2a32 f(ti−1 . con lo que el sistema reducido de ecuaciones nos queda ⎧ ⎪ b1 + b2 + b3 = 1. Entonces la quinta y sexta ecuaciones se simpliﬁcan a b2 c2 + b3 c2 = 2 3 1 3 que es la cuarta. ⎪ 2 3 ⎪ 3 ⎩ b3 a32 c2 = 1 . obtenemos las ecuaciones ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ b1 + b2 + b3 = 1. y comparando con el desarrollo de Taylor de orden 2 (2. 21 3 b3 a32 c2 = 1 . yi−1 )f(ti−1 . yi−1 ) + a21 f(ti−1 . 6 que nos dan los métodos de Runge—Kutta de orden 3.10).9). llegamos a c3 = a31 + a32 . que es una familia biparamétrica de métodos numéricos. 2 3 3 b2 c2 a21 + b3 c3 (a31 + a32 ) = 1 . yi−1 ) + 2c3 (a31 + a32 ) f(ti−1 . yi−1 ) f(ti−1 . yi−1 ) h2 ∂y ∂t ∂y +O(h2 ) Introduciendo los desarrollos de Taylor de G2 y G3 en la ingualdad (2. ⎪ b2 c2 + b3 c3 = 1 . b2 c2 + b3 c3 = 1 . 2 b2 c2 + b3 c2 = 1 . ⎪ ⎪ ⎪ c2 = a21 . y la última ecuación es la antepenúltima.
dado que el error local de los métodos de Runge—Kutta de tres etapas es de orden 4 [O(h4 )]. = 6b3 c2 6 1 − c22 3 y de la última De c3 = a31 + a32 tenemos que a31 = c3 − a32 = c3 − y de b1 + b2 + b3 = 1 concluimos b1 Así.6. ¯ 2 2 ¯ ¯ c2 c3 ¯ a32 = 1 c2 c3 (c3 − c2 ) ¡ ¢ . 2. aumentando en un orden los desarrollos de Taylor de las funciones implicadas. Existen métodos que se obtienen cuándo alguna de estas cantidades son nulas. tenemos que si la tabla de Butcher de un método de cuatro etapas es 0 c2 c3 c4 a21 a31 a41 b1 a32 a42 b2 43 a43 b3 b4 . Como vemos. 6 1 − c22 3 0 2 9 3 4 1 3 4 9 y 0 1 2 1 2 1 −1 1 6 2 2 3 1 6 Estas soluciones son válidas siempre que c2 y c3 sean no nulos y distintos.Métodos de un paso ¯ ¯ ¯ c2 1 ¯ 2 ¯ 2 1 ¯ ¡ ¢ ¯ ¯ c2 c2 1 − c22 3 ¯ 3 b3 = ¯ ¯ c2 c3 ¯ = c2 c3 (c3 − c2 ) . tenemos que el error global de los métodos de tres etapas es de orden 3.2 El método de cuatro etapas Procediendo como en los casos de dos y tres etapas. obtenemos los métodos de Runge—Kutta de orden tres 0 1 2 3 4 1 2 ¡ ¢ ¡ ¢ c3 c23 − 1 c2 1 − c22 3 3 = 1 − b2 − b3 = 1 − − c2 c3 (c3 − c2 ) c2 c3 (c3 − c2 ) ¡ ¢ ¡ ¢ c2 c3 (c3 − c2 ) − c3 c23 − 1 − c2 1 − c22 3 3 = c2 c3 (c3 − c2 ) ¡ ¢ c2 c3 + 1 (c3 − c2 ) − 1 (c2 − c2 ) 2 3 2 3 = c2 c3 (c3 − c2 ) c2 c3 (c3 − c2 ) ¡ ¢ . y que se obtienen de igual manera.
(2.12) obtenemos a43 = (1 − c2 )(2c2 − 1)(1 − c3 ) . de la cual se tiene para j = 4 que b4 (1 − c4 ) = 0. 12c2 (1 − c2 )(1 − c3 ) Tomando ahora j = 3 en (2. 4. ⎪ 8 ⎪ ⎪ ⎪ b3 a32 c2 + b4 (a42 c2 + a43 c2 ) = 1 . de donde c4 = 1 ya que b4 6= 0.12) con j = 2 obtenemos a32 = y a42 = (1 − c2 )[2(1 − c3 )(1 − 2c3 ) − (c2 − c3 )] . j = 2. 12c3 (1 − c3 )(c3 − c2 ) 6c2 c3 − 4(c3 + c2 ) + 3 . ⎪ 2 2 3 ⎪ 12 ⎩ 1 b4 a43 a32 c2 = 24 . ⎪ 2 ⎪ ⎪ ⎪ b2 c2 + b3 c2 + b4 c2 = 1 . bi aij = bj (1 − cj ). 3. j = 2. b3 y b4 . 6 1 3 3 3 ⎪ ⎪ b2 c2 + b3 c3 + b4 c4 = 4 . Butcher en 1963 dio la simpliﬁcación 4 X i=1 ⎧ ⎪ b1 + b2 + b3 + b4 = 1. 2c2 (2c2 − 1) Finalmente.12) donde aij = 0 si i ≥ j. como sabemos. 4. 3 y 5 son lineales en b2 . ⎪ 2 3 4 3 ⎪ ⎨ b3 a32 c2 + b4 (a42 c2 + a43 c3 ) = 1 . la primera de las cuales. 2c2 (c2 − c3 )[6c2 c3 − 4(c2 + c3 ) + 3] 44 . 3.Métodos de un paso entonces han de satisfacerse la simpliﬁcación cj = y las condiciones j−1 X k=1 ajk . ⎪ ⎪ b3 c3 a32 c2 + b4 c4 (a42 c2 + a43 c3 ) = 1 . de la última ecuación y (2. Como las ecuaciones 2. viene de la condición de consistencia de los métodos de Runge—Kutta en general. tenemos que b2 = b3 = y b4 = 1 − 2c3 . c3 (c2 − c3 )(6c2 c3 − 4(c2 + c3 ) + 3) c3 (c2 − c3 ) . ⎪ ⎪ ⎪ b2 c2 + b3 c3 + b4 c4 = 1 . 12c2 (1 − c2 )(c2 − c3 ) 1 − 2c2 .
1/2.3 Métodos de más etapas En general. De hecho. Cuando alguna de estas condiciones no se satisfacen. a partir de 5 etapas no aumenta el orden global del método. 1}. 0 1 2 1 2 −1 3 1 1 8 1 −1 3 8 1 3 8 1 8 1 2 1 0 0 1 6 1 2 0 1 3 1 1 3 1 6 que es el método debido a Kutta de 1905 y que usualmente se conoce como el método de Runge— Kutta. B = log C + p log(tf − t0 ). / / c3 ∈ {0. 2. se tiene la siguiente tabla para los órdenes en función de las etapas Etapas 1 2 3 4 5 6 7 8 Orden 1 2 3 4 4 5 6 6 es decir. Para ﬁjar ideas. existen no obstante métodos de Runge—Kutta con estos coeﬁcientes. donde p es el orden. ¿Por qué entonces se construyen métodos de más de 5 etapas? La razón es meramente computacional. supongamos que el error global de un método es e = Chp . Tomando logaritmos tenemos que log e = log C + p log h.Métodos de un paso Como vemos. y suponiendo que el tamaño del paso h = (tf − t0 )/n se tiene log e = log C + p log tf − t0 n = log C + p log(tf − t0 ) − p log n = B − p log n. 45 siendo . Estos métodos son de orden 4 (orden 5 tiene el error local de truncamiento) y representan una familia de cinco parámetros de métodos cuyos ejemplos son 0 1 3 2 3 1 3 1 y sobre todo. se conoce que un método de Runge—Kutta de m etapas tiene a lo sumo orden m. 1} y c2 6= c3 .6. estas soluciones dependientes de dos parámetros son válidas siempre que c2 ∈ {0. que cumple c2 = c3 . ya que los errores debidos al redondeo aumentan con el número de pasos.
Métodos de un paso Si asumimos cierta la identidad log e = B − p log n y escirbimos otra análoga para un método de orden q log e = B 0 − q log n. 46 . tenemos dos rectas que se cortarán en un punto. que marcarán la eﬁciencia de cada método y determinarán cuándo debe usarse cada uno.
. bj y p son parámetros elegidos de acuerdo con unas condiciones de convergencia y estabilidad detrminadas. . Convergencia y estabilidad de los métodos multipaso. t0 ..1 Introducción Partimos de la ecuación diferencial con condiciones iniciales ½ 0 y = f(t.. Métodos de Adams—Moulton. y(t0 ) = y0 . Si p ≥ 1. y−j ). Métodos de Adams—Bashforth. y).1) donde aj .. y sean tf = t1 > t0 .. 1.Capítulo 3 Métodos multipaso Sumario. 3. p. Métodos multipaso generales. Si b−1 = 0. que veriﬁca tener unicidad de solución.. entonces para j = 1.2 Métodos multipaso lineales Sea y(t. y0 ). y0 ) la solución del problema de condiciones iniciales anterior. es decir. Métodos predictor—corrector. Hasta ahora hemos visto métodos numéricos de un paso. e y−j ≈ y(t0 − jh. la sucesión que aproxima la solución yi se genera de forma recursiva a partir de los términos inmediatamente anteriores. j = 0. Como veremos. en los métodos multipaso esta sucesión se construye a partir de una ecuación en diferencias con orden mayor que uno. (3. Sea y1 la aproximación de y(t1 . t0 . 3. p los valores y−j han de ser previamente estimados con un método de orden uno.. y h = t1 − t0 . el método se dirá 47 . Veamos a continuación cómo construir tales aproximaciones. y0 ) dada por la expresión y1 = p X j=0 aj y−j + h j=−1 p X bj f(t0 − jh. probablemente algún método de Runge—Kutta. t0 . esto es yi−1 .
y0 ) (k − 1)! j=−1 k! j=0 k! 48 Ã j=−1 j=−1 X p jbj + 1X 2 1 j aj − 2 j=0 2 p ! j=0 . y0 ) − y(ti + h. Veamos a continuación cómo obtener los parámetros del método. . t0 . t0 .. k! k=0 X ∞ X k=0 1 k) y (ti .. consideramos el error local de truncamiento p p X X ti+1 = aj y(ti − jh. (3. aj p ∞ X 1 +h bj yk+1) (ti .2) aj y(ti − jh.. y0 )(−1)k (−ih)k k! y reagrupando en distintas potencias de h concluimos que Ã p ! Ã p ! p X X X aj − 1 + hy0 (ti . t0 . j bj + j aj − +. y−j ).. y0 ) + h bj f(ti − jh. n} construimos yi+1 = p X j=0 aj yi−j + h j=−1 junto con las condiciones yj .. 1. y0 ) − Ã ! p p 1 X 2 1X 3 1 j bj − j aj − +h3 y3) (ti . y0 ) j=0 +h2 y00 (ti . t0 . t0 . . y0 )hk . y(ti − ih. t0 . y0 ) bj − iaj − 1 ti+1 = y(ti .. y0 ) = X j=0 j=0 p j=−1 p p X bj f(ti − jh. y0 )(−1)k (−ih)k k! j=−1 k=0 ∞ X 1 − yk) (ti . y0 ) − y(ti + h. p.. j = 0.1). y0 ). y0 ) = X j=0 j=0 p j=−1 X bj y0 (ti − jh. t0 .Métodos multipaso explícito ya que y1 se obtiene directamente como p p X X y1 = aj y−j + h bj f(t0 − jh. t0 . y0 )) − y(ti + h. tenemos p p X X ti+1 = aj y(ti − jh. t0 . h = (tf − t0 )/n y ti = t0 + ih. para cada i ∈ {1. t0 . t0 . j=0 j=0 Sin embargo. y0 ) + h j=−1 y tomando el desarrollo en serie de Taylor en h = 0.. En general. t0 .. y0 ) 2 j=−1 6 j=0 6 ! Ã p p (−1)k−1 X k−1 (−1)k X k 1 k k) + . t0 . Para ello. Como vemos se trata de una ecuación en diferencias de orden p que da lugar a la aproximación de la solución. t0 . yi−j ). presumiblemente haciendo uso de algún método numérico para ello. si b−1 6= 0 entonces el método se dice implícito porque hay que calcular y1 resolviendo la ecuación (3. + h y (ti .. t0 . y0 ) + h bj y0 (ti − jh. t0 . t0 . si tenemos n pasos.
a1 = 5. yi ) + 2hf(ti − h.. yi−1 ). q.. satisfaciendo el sistema de ecuaciones ⎧ ⎪ a0 + a1 = 1. k = 1. ⎪ ⎨ −a1 + b−1 + b0 + b1 = 1. ⎪ ⎪ b0 = 2 − 3λ. Será entonces un método de la forma yi+1 = a0 yi + a1 yi−1 + hb0 f(ti . yi−1 ). (−1) k j=−1 Veamos algunos ejemplos concretos de métodos multipaso. yi+1 = −4yi + 5yi−1 + 4hf(ti . de donde b1 = 2. podemos escribir de forma más compacta Ã p ! Ã ! p ∞ X X 1 X ti+1 = y(ti . que es un sistema compatible indeterminado cuya solución es ⎧ ⎪ a0 = 1 − λ. Supongamos en primer lugar que el método es explícito. 2. ⎪ ⎩ −a1 + 3b−1 + 3b1 = 1. λ ∈ R. esto es. esto es b−1 = 0. yi+1 ) + hb0 f(ti . ⎪ a1 − 2b1 = 1. el sistema de ecuaciones anterior es de la forma ⎧ ⎪ a0 + a1 = 1.. 3. ⎨ 5 b−1 = 12 − λ.3 Primeros ejemplos Vamos a construir un método multipaso de 2 pasos y de orden 3. ⎪ ⎩ −a1 + 3b1 = 1. ⎪ a1 + 2b−1 − 2b1 = 1. ⎪ ⎨ −a1 + b0 + b1 = 1. k! j=0 j=−1 k=1 con lo que el error local de truncamiento será de orden O(hq+1 ) si se satisfacen las ecuaciones p X j=0 aj = 1.Métodos multipaso y si ﬁjamos a−1 = 0. y0 ) hk yk) (ti . yi ) + hb1 f(ti − h. t0 . de donde Si elegimos un método implícito. yi ) + hb1 f(ti − h. yi−1 ). p X j k−1 (jaj − kbj ) = 1. . ⎪ ⎪ ⎪ a1 = λ. t0 . y0 ) (−1)k aj − 1 + j k−1 (jaj − kbj ) − 1 . 49 . yi+1 = a0 yi + a1 yi−1 + hb−1 f(ti + h. b0 = 4 y a0 = −4. ⎪ 3 ⎪ ⎩ 1 b1 = − 12 + 5λ.
Sin embargo. el sistema dado por no tiene solución. y tomando las últimas 4 ecuaciones. y(0) = 1. yi−1 ). cosa que no ocurre con el método explícito. el ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ determinante de la matriz ampliada sería ¯ 1 2 −2 1 ¯ ¯ −1 3 3 1 ¯ ¯ = −8. la ecuación adicional a1 − 4b1 = 1. tenemos el método yi+1 = yi + h 5 2 1 f(ti + h. Por ejemplo. yi+1 ) + h f(ti . −a1 + 3b1 = 1. si consideramos el problema ½ 0 y = y. yi ) − h f(ti − h. ya que debería cumplir. 12 3 12 El método implícito tiene la desventaja de que el valor de yi+1 ha de obtenerse a partir de la solución aproximada de una ecuación algebraica. y dado que el determinante ¯ ¯ ¯ 1 −2 1 ¯ ¯ ¯ ¯ −1 3 1 ¯ = 4. Este último no permite obtener aproximaciones de tanto orden como el implícito. 1 4 −4 1 ¯ ¯ −1 5 5 1 ¯ por lo que el sistema será incompatible y no habrá solución del mismo. Así pues. Ahora bien. Si intentáramos conseguir orden 5. el método de dos pasos explícito no puede tener orden 4. y ahora el determinate ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ 1 1 0 0 0 0 −1 1 1 1 0 1 2 0 −2 0 −1 3 0 3 0 1 4 0 −4 ¯ ¯ ¯ ¯ ¯ ¯ = −12. tendríamos que añadir la ecuación −a1 + 5b−1 + 5b1 = 1. ¯ ¯ ¯ ¯ por lo que el sistema lineal anterior es compatible determinado y tiene un solución única que será el único método de orden 4 y paso 2. los métodos de dos pasos pueden tener a lo sumo orden 3 si es explícito y orden 4 si es implícito. si planteamos el mismo sistema para el método implícito.Métodos multipaso por lo que si λ = 0. además de las ecuaciones anteriores. ⎩ a1 − 4b1 = 1. ¯ ¯ ¯ 1 −4 1 ¯ ⎧ ⎨ a1 − 2b1 = 1. 50 . hemos de añadir la ecuación a1 + 4b−1 − 4b1 = 1.
Integrando respecto de la variable independiente Z ti+1 Z ti+1 0 y(ti+1 . y0 )dt = f(t. ti+1 = y(ti . . 51 . y(ti − jh. t0 ..Métodos multipaso y queremos aproximar y(1. i.. y0 ) = − E(t)dt. 0. n. y(t.. j = i − p. t0 . Partamos de la ecuación diferencial y buscamos aproximar y(tf . y0 ) − y(ti . yi+1 = yi + ti donde es decir. yi−j ). t0 . el polinomio interpolador sobre los datos anteriores que serán conocidos. dependerá del error que se comete al calcular el polinomio interpolador. y que cumplirá la condición Q(ti − jh) = f(ti − jh. t0 . . y0 = 1. y0 )) = Q(t) + E(t). ti ti Si queremos cosntruir un método de p pasos. el error global crece de forma dramática. y1 . pero antes veamos cómo deducir los métodos de multipaso a partir de la integración numérica. donde y1 lo hemos elegido a partir de un método de Runge—Kutta de orden 3 para un tamaño de paso preﬁjado h..4 Métodos de multipaso deducidos a partir de la integración numérica ½ y0 = f(t. y0 ))dt. vemos que la aproximación empieza a oscilar. sustituimos f(t. y(t. 1). y0 ). t0 .. que denotaremos Q(t). t0 . y0 ) = y (t. y). t0 . con lo cual. . j = i − p. El error local de truncamiento será en este caso Z ti+1 Z ti+1 Q(t)dt − y(ti+1 . a pesar de tener una aproximación local de orden 3. y que será clave a la hora de determinar el error local cometido en la aproximación. El método numérico que hemos de construir tendrá la forma Z ti+1 q(t)dt.. i. i = 0. Vamos a ver a continuación cómo se construyen los métodos de multipaso explícitos e implícitos. Este polinomio de interpolación veriﬁca que f(t. t0 . y0 ) + ti ti q(ti − jh) = f(ti − jh.. y aplicamos el método ½ yi+1 = −4yi + 5yi−1 + 4hyi + 2hyi−1 . y0 )). 1. es decir. t0 ... y(t0 ) = y0 . y(t. 3. Veremos qué ocurre con este método. Para ello ﬁjamos h = (tf − t0 )/n y ti = t0 + hi. esto es tf = 1. donde E(t) es el error cometido al aproximar la función por el polinomio interpolador. t0 . y0 )) por un polinomio de interpolación de grado p.
= j j=0 µ s+p−1 p ¶ ∇p f(ti . se deﬁnen ∇1 xi = ∇xi = xi − xi−1 . yi ) donde µ m k ¶ = m(m − 1). t0 .(m − k + 1) k! µ m 0 ¶ para m ∈ R y k ∈ N.Métodos multipaso 3. + 2 p Xµ s + j − 1 ¶ ∇j f(ti . Si t = ti + sh. el error del polinomio de interpolación E(t) = Como s = t−xi .. y ∇3 xi = ∇2 xi − ∇2 xi−1 = (xi − 2xi−1 + xi−2 ) − (xi−1 − 2xi−2 + xi−3 ) = xi − 3xi−1 + 3xi−2 − xi−3 .. Por ejemplo ∇2 xi = ∇xi − ∇xi−1 = (xi − xi−1 ) − (xi−1 − xi−2 ) = xi − 2xi−1 + xi−2 . y para k > 1. y0 )) + O(hp+2 ). estableceremos que ∇0 xi = xi . Además. y(ti .4. h µ s+p p+1 ¶ hp+1 f p+1) (ti .. Dada la sucesión xi . el polinomio de interpolación tiene la forma q(t) = f(ti . yi ) + s∇f(ti . 0 . yi ).. Por convenio. yi ) + . y = 1. reescribimos yi+1 = yi + Z ti+1 q(t)dt = yi + h 52 ti Z 1 q(xi + sh)ds.1 Métodos de Adams—Bashforth Para cosntruir estos métodos utilizamos la forma de Newton del polinomio interpolador dada por la diferencias ﬁnitas. yi ) + µ s+1 2 ¶ ∇ f(ti . ∇k xi = ∇k−1 xi − ∇k−1 xi−1 .
. siendo el método de la forma µ ¶ 1 5 2 yi+1 = yi + h f(ti . pero ahora tomando el polinomio interpolador en t = ti+1 + sh. 12 Como vemos. Así. j = 1. con lo que ¶ ¶ µ µ s+1 s+p−1 2 ∗ ∇ f(ti . p. t0 .Métodos multipaso El método de Adams—Bashforth se construye a partir del desarrollo Z 1 yi+1 = yi + h q(xi + sh)ds 0 ÃZ ! ¶ p 1 Xµ s+j−1 ∇j f(ti . yi+1 ) q (t) = f(ti+1 . + ∇p f(ti+1 . ¶ Z 1µ s+j−1 ds. yi−2 )] 12 1 = yi + h (23f(ti .2 Métodos de Adams—Moulton Los métodos de Adams—Moulton se obtienen de igual manera que los explícitos de Adams—Basfhforth. yi ) + ∇ f(ti . yi ) 2 12 µ 1 = yi + h f(ti . necesitaremos que p = 2. los primeros coeﬁcientes son γ1 1 2 X j=0 j=0 p γ j ∇j f(ti . 2. yi ) − 16f(ti−1 . yi ) j 0 = yi + h donde los coeﬁcientes γ 0 = 1. yi−1 ) + 5f(ti−2 . yi ) − 2f(ti−1 . 3. yi ) + f(ti−1 . yi+1 ) + 2 p 53 . y0 ) + O(hp+3 ). por lo que será convergente. yi ) + ∇f(ti ...4. yi )ds = yi + h j 0 j=0 Ã p ¶ ! Z 1µ X s+j−1 ds = yi + h ∇j f(ti . yi ) + [f(ti . yi−1 )] 2 ¶ 5 + [f(ti . si buscamos un método de error global O(h3 ). γj = j 0 se calculan de forma directa. De hecho. yi+1 ) + s∇f(ti+1 . . yi ) + .. yi−1 ) + f(ti−2 . γ2 5 12 γ3 3 8 γ4 251 720 95 288 γ5 γ6 19087 60480 Además.. por lo que es orden O(hp+2 ). yi ). al método de Adams—Bashforth le podemos aplicar el teorema de convergencia global. el error local de truncamiento es de la forma ti+1 = −hp+2 γ p+1 yp+2 (ti . yi−2 )) . al contrario de lo que ocurría con el primer ejemplo que estudiamos.
como s = t−xi+1 . yi+1 )ds = yi + h j −1 j=0 Ã p ¶ ! Z 0µ X s+j−1 ds ∇j f(ti+1 . yi+1 ) − ∇ f(ti+1 . yi+1 ). yi ) − f(ti−1 .. yi+1−j ).. 2. j β ∗ f(ti+1−j . yi−1 )) .. j = 1. = yi + 12 A la hora de aproximar yi+1 . yi+1 ). = j j=0 con lo que. y0 ) + O(hp+3 ). yi+1 ) 2 12 h (5f(ti+1 . j 54 . j γ ∗ = 1. . yi+1 ) = yi + h j −1 0 Z 0 q(xi + sh)ds ÃZ −1 = yi + h donde X j=0 j=0 p γ ∗ ∇j f(ti+1 .Métodos multipaso p Xµ s + j − 1 ¶ ∇j f(ti+1 . yi+1 ) + 8f(ti . 0 ¶ Z 0µ s+j−1 ∗ ds. h yi+1 Desarrollando. p. γj = j −1 y el error local de truncamiento es ti+1 = −hp+2 γ ∗ yp+2 (ti . yi+1 ). t0 . yi+1 ) − ∇f(ti+1 . obtenemos reescribimos Z ti+1 Z ∗ = yi + q (t)dt = yi + h ti 0 q∗ (xi + sh)ds. −1 yi+1 = yi + h ! p Xµ s + j − 1 ¶ ∇j f(ti+1 . démonos cuenta de que si deﬁnimos G(yi+1 ) = yi + h = yi + h p X X j=0 j=0 p γ ∗ ∇j f(ti+1 . p+1 Los primeros coeﬁcientes son en este caso γ∗ 1 −1 2 γ∗ 2 1 − 12 γ∗ 3 1 − 24 γ∗ 4 19 − 720 γ∗ 5 3 − 160 863 − 60480 γ∗ 6 Tomando p = 2 obtenemos el método de dos pasos y orden tres µ ¶ 1 1 2 yi+1 = yi + h f(ti+1 .
0 3.. y0 ) la solución del problema de condiciones iniciales ½ 0 y = f(t.. Para ello. − ap yi−p = 0. t0 . z)|| 0 ≤ h|β ∗ |L||y − z||. (3.. y(t0 ) = y0 . y(0) = y0 . sabemos que p(1) = 1 − p X j=0 p X j=0 aj yi−j . − ap . aj = 0. Como sabemos. y h = t1 −t0 . Vamos a ver qué propiedad adicional hemos de añadir a la convergencia local para que el método sea globalmente convergente. k ≤ p+1.Métodos multipaso por lo que ||G(y) − G(z)|| = h|β ∗ |||f(ti+1 . que es una ecuación en diferencias lineal cuyo polinomio característico es p(λ) = λp+1 − a0 λp − . Así. y).. y) − f(ti+1 . y sean tf > t0 . n Sea yi+1 la aproximación de y(tt+1 . yi−j ). 0 dado que f es Lispchitziana en la variable y. con lo que el método multipaso queda como yi+1 = o equivalentemente yi+1 − a0 yi − a1 yi−1 − . aplicamos el método al problema ½ 0 y = 0. y0 ) dada por el método multipaso yi+1 = p X j=0 aj yi−j + h p X j=0 bj f(ti − jh. Por la condición de convergencia local. 55 . para que se pueda aplicar el Teorema del punto ﬁjo a G hemos de elegir tamaños de paso h suﬁcientemente pequeños para que h|β ∗ |L < 1.3) donde los coeﬁcientes se han escogido de manera que el error local de truncamiento es O(hk ). la convergencia local no implica necesariamente la convergencia global.5 Estabilidad de los métodos multipaso lineales Sea y(t. t0 .
0 a0 1 0 .. .. En este caso. ⎪ p ⎪ zi+1 = zp+1 . ⎟ = ⎜ . que tiene por raíces 1 y −5.. Dicha ecuación.. Para conseguir un método convergente.. zi+1 . . yi−1 )... p. ⎪ 2 ⎪ z = z3 . ⎜ p ⎟ ⎜ ⎝ zi+1 ⎠ ⎝ 0 a0 zp+1 i+1 tiene radio espectral ρ(A) = max {|μ| : μ valor propio de A} = max {|μ| : p(μ) = 0} ≤ 1... zi = 0. vemos que p(λ) = λ2 + 4λ − 5. la estabilidad también implica la convergencia.Métodos multipaso por lo que 1 es solución particular de la ecuación en diferencias. quedarnos con las ecuaciones ⎧ ⎨ a0 + a1 = 1.3) puede escribirse como ⎧ 1 2 ⎪ zi+1 = zi . . y si éste es uno. 0 ap−1 0 0 ... ⎪ i+1 ⎨ i . Si b−1 6= 0. ⎩ a1 − 2b1 = 1. será entonces estable si las restantes raices de p(λ) tienen módulo menor o igual que uno..... por lo que el método no es estable. 1 ap ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎟ ⎟ ⎟.. y de ahí su divergencia.. entonces si b−1 = 0.. 0 a2 lo escribimos como ⎛ ⎞ 0 z1 i ⎜ 0 z2 ⎟ i ⎜ ⎟ ⎟+ h⎜ .. aunque la demostración sale fuera de los contenidos del curso....... Si consideramos el método del ejemplo inicial. 0 0 0 1 . ⎟ · ⎜ ⎟ ⎜ 0 0 .. zi+1 = j=0 aj zp−j+1 + h p i j=−1 bj f(ti − jh.. 0 1 ⎠ ⎝ a1 a2 . 1. ⎟ ⎜ p ⎠ ⎝ 0 zi Pp p−j+1 p+1 ) zi j=−1 bj f(ti − jh.. j i ⎛ 1 ⎞ ⎛ zi+1 0 ⎜ z2 ⎟ ⎜ 0 ⎜ i+1 ⎟ ⎜ ⎜ ... . 0 a1 0 1 . el método (3.. ..... 56 . .. 0 0 ⎟ ⎜ ⎟ ⎜ .... se trata de una raíz simple. . el método multipaso veriﬁca las condiciones del teorema de convergencia global. ⎪ i ⎪ p+1 Pp P ⎩ p−j+1 ).. Matricialmente ⎞ ⎛ 1 0 . 0 0 . ⎟ ⎠ Si la matriz donde zp−j+1 = yi−j .... . yi ) + 2hf(ti − h.. esto es. por lo que si el método es estable de orden local de truncamiento O(hk+1 ). .. hemos de imponer que el orden local de truncamiento sea una unidad menor.... ..... ap−1 ap ⎜ ⎜ A =⎜ ⎜ ⎝ ⎛ 0 0 .. entonces será convergente de orden O(hk+1 ). yi+1 = −4yi + 5yi−1 + 4hf(ti .. −a1 + b0 + b1 = 1.
q(ti+1 )). t0 . j = t −t 0. ⎪ ⎩ b1 = μ. q(s) = j j=0 cuyo error es E(t) = y suponemos que se veriﬁca la ecuación diferencial para dicho polinomio en ti+1 . conocido como la regla del punto medio. ⎪ ⎨ a1 = 1 + 2μ. y0 )) + O(hp+2 ). de la solución aproximada dada por p Xµ s + j − 1 ¶ ∇j yi+1 .. yi+1 = yi−1 + 2hf(ti .. Por una parte. yi ) + μhf(ti−1 . −1 ≤ μ ≤ 0.6 Fórmulas BDF Las fórmulas BDF se basan en la utilización del polinomio interpolador en los puntos yi+1−j . obtenemos el método p(λ) = λ2 + 2μλ − (1 + 2μ). . ⎪ b0 = 2 + μ. . y por otra.. −1 ≤ 1 + 2μ ≤ 1. 3. q(ti+1 ) = yi+1 . y(ti . esto es. yi−1 ). 1. p.Métodos multipaso que nos da el conjunto uniparamétrico de soluciones ⎧ ⎪ a0 = −2μ. que da lugar a la familia de métodos yi+1 = −2μyi + (1 + 2μ)yi−1 + (2 + μ)hf(ti . con polinomio característico que para ser estable de cumplir que o equivalentemente y eligiendo μ = 0. yi ). dado que t = ti+1 + sh ! Ã p µ X s+j−1 ¶ d q0 (s) = ∇j yi+1 j ds j=0 µ ¶ p X s+j−1 −1 d ∇j yi+1 . Es obviamente estable ya que p(λ) = λ2 − 1. μ ∈ R. h = j dt j=0 57 µ s+p p+1 ¶ hp+1 f p+1) (ti . h = f n 0 . cuyas raices son ±1. q0 (0) = f(ti+1 .
Por ejemplo. el tamaño de paso h debe ser elegido para que el método iterativo sea convergente con constante de Lipschitz menos que uno... yi+1 ) = 2 X1 j=1 p X1 j=1 j ∇j yi+1 = hf(ti+1 . hf(ti+1 . y como corrector. 58 (3. λ ∈ R. Por ejemplo. consideramos de entre la familia yi+1 = a0 yi + a1 yi−1 + hb−1 f(ti + h. 1 ∇j yi+1 = ∇yi+1 + ∇2 yi+1 j 2 1 (yi+1 − 2yi + yi−1 ) 2 3 1 yi+1 − 2yi + yi−1 . de igual manera que pasaba con los métodos de Adams implícitos. puede probarse que es de orden O(hp+1 ). j! ½ (j−1)! j! se tiene que 0 gj (0) = = 1 j 0 si j > 0. ⎪ ⎪ b0 = 2 − 3λ.Métodos multipaso Como µ ¶ s+j−1 j = (s + j − 1)(s + j − 2). ⎨ 5 b−1 = 12 − λ. yi−1 ). = 2 2 = yi+1 − yi + y como vemos es de dos pasos. 3. Como estos métodos son implícitos. yi+1 ) + hb0 f(ti . de donde que se conocen como las fórmulas BDF. ⎪ ⎪ ⎪ a1 = λ. yi+1 ). si j = 0. El método explícito se usa para obtener la condición inicial con la que obtener el mediante un método iterativo.7 Metodos predictor—corrector Los métodos de predictor corrector se basan en utilizar alternativamente métodos multipaso explícitos e implícitos de un mismo orden para aproximar la solución. si p = 2. yi ) + hb1 f(ti − h. En cuanto al error de truncamiento local. una mejor aproximación con el método implícito.4) . ⎪ 3 ⎪ ⎩ 1 b1 = − 12 + 5λ.(s + 1)s = gj (s). dados por el sistema ⎧ ⎪ a0 = 1 − λ. consideramos como predictor el método explícito yi+1 = −4yi + 5yi−1 + 4hyi + 2hyi−1 .
en general utilizaríamos un método de Runge—Kutta.5) sí es estable y por tanto convergente.4). tomando como punto ∗ inicial para hacer las iteraciones el punto y2 . • Cuando el valor obtenido de y2 sea aceptable. 59 . ya que el coste computacional de cálculo de f es menor en este caso. 12 3 12 (3. ∗ • Predecimos el valor de y2 por y2 con el método (3.5) será convergente. • Mejoramos la aproximación anterior calculando y2 con el método (3. por lo que el método con λ = 0. las distintas evaluaciones de f en etapas anteriores se pueden utilizar en varios pasos. utilizamos un método de Runge—Kutta de orden 3 para estimar y1 . = 2 2 1−λ± que nos da 1 y −λ como soluciones. yi+1 ) + h f(ti . Como el error local es de orden 4. esta aproximación será de este orden.Métodos multipaso Para obtener un método convergente. Un esquema para aplicar estos métodos sería el siguiente: • Como el método es de orden 2. por lo que no cada vez que se aplica el método multipaso hay que evaluar la función f en todos los pasos. Hemos de destacar que el método (3. yi ) − h f(ti − h. 3.8 Multipaso o Runge—Kutta La elección del método numérico utilizado para obtener la solución aproximada depende en gran medida del coste de computación de la función f. yi−1 ).5). que son p (1 − λ)2 + 4λ t = p2 1 − λ ± (1 + λ)2 1 − λ ± (1 + λ) = . y así sucesivamente. Si este coste es bajo. calculamos las raices de p(t) = t2 + (λ − 1)t − λ. volvemos a aplicar los dos puntos anteriores para obtener y3 . dado por yi+1 = yi + h 5 2 1 f(ti + h. Además. pero cuando este es alto utilizamos uno de multipaso.
Métodos multipaso 60 .
Parte II Bloque de práctica 61 .
por lo que esta práctica trata de explicar como pedirle a Mathematica que haga aquello que nosotros deseamos. gtraﬁca de funciones. Con Mathematica ocurre lo mismo. Además. Mathematica es un programa que permite hacer cálculos matemáticos complicados con gran rapidez. Por ejemplo. Funciones. Pretendemos con esta práctica introducir al alumno en el manejo de este potente programa que puede servirle de utilidad en futuros cálculos que deba realizar. Sólo conociendo teóricamente las operaciones que Mathematica realiza nos percataremos de su utilidad. Generalidades de Mathematica. Abordaremos en esta práctica una iniciación a Mathematica partiendo desde cero. hemos de pulsar las teclas shift + enter para obtener el resultado. 4. debemos hacer hincapié en el hecho de que es necesario por parte del alumno un conocimiento matemático teórico de todas las funciones y sentencias que vamos a usar. sino que también hace derivadas. cálculo de primitivas. intentando poner de maniﬁesto su utilidad a la hora de trabajar con expresiones matemáticas complicadas. es como una calculadora gigante a la que no sólo podemos pedirle que haga cálculos numéricos. aunque una calculadora multiplica números con suma facilidad. bien sean éstas numéricas o simbólicas. Por ejemplo. aparece una pantalla blanca vacía. representación gráﬁca de curvas y superﬁcies. derivadas. etcétera. sólo nos percatamos de su potencia en cuanto conocemos dicha operación y somos capaces de realizarla de un modo mucho más lento. representación . el programa puede utilizarse para corregir los problemas propuestos al alumno en la clase de problemas. A pesar de la utilidad del programa. permitiendo hacer operaciones en poco coste de tiempo y con bastante facilidad.1 Preliminares Cuando se arranca Mathematica.Capítulo 4 Introducción a Mathematica Sumario. En ella podemos escribir aquellas operaciones que queremos que realice. supongamos que queremos hacer la 63 . Esencialmente vamos a aprender a utilizar el programa. Una vez tecleada la operación. Para entendernos.
operará estos números ofreciendo siempre su valor exacto. Teclearemos entonces 2+2 en la pantalla.Introducción a Mathematica operación 2 + 2. Además de la suma se pueden realizar las siguientes operaciones algebraicas como si se tratara de una calculadora: x+y suma de números x−y resta de números x/y división de números x y x ∗ y producto de números xˆy potencia xy Cuando Mathematica realiza alguna de las siguientes operaciones. n].619047619 64 .619048. se tiene In[2] := 1/3 + 2/7 13 . Para ello se tienen las sentencias x//N N[x] N[x. Las primeras escriben el número x con seis cifras signiﬁcativas. como veremos posteriormente. si escribimos In[3] := 1/3 + 2/7 //N Out[3] = 0. es decir. Podrá aparecer únicamente un nú mero de entrada.0 del programa y posteriores esta última sentencia no siempre funciona del modo deseado). 10] Out[4] = 0. Si por el contrario escribimos In[4] := N[1/3 + 2/7. obtendremos el resultado con 6 cifras signiﬁcativas. a veces nos es más útil tener el valor de este número expresado con cifras decimales. Out[2] = 21 Sin embargo. Al ir explicando las diferentes operaciones que Mathematica realiza. iremos escribiéndolas en la forma en que el programa lo escribe en la pantalla de ordenador. por ejemplo 1/3 + 2/7. A continuación pulsamos mayúsculas + enter o la tecla intro en el teclado numérico y a continuación aparecerá en pantalla In[1] := 2 + 2 Out[1] = 4. mientras que la segunda escribe dicho número con un número n de cifras signiﬁcativas que nosotros preﬁjamos (en la versión 4. Por ejemplo. Todas las operaciones realizadas por el programa cuando se pulsan las teclas mayúsculas + enter tienen asignadas un número de entrada marcado por In[·] y el mismo número de salida cuando se realiza la operación marcado por Out[·].
bien sea una operación especíﬁca del programa. Hemos de señalar que Mathematica sigue el orden conocido de preferencia sobre las operaciones. como por ejemplo matrices o vectores. • Las llaves {·} se utilizan para asignar valores numéricos a las variables. A grosso modo podemos indicar las siguientes generalidades: • Los paréntesis se usan en las operaciones algebraicas para indicar la preferencia a la hora de hacer las operaciones. En caso de las operaciones numéricas también tendremos una valor numérico aproximado con seis cifras signiﬁcativas si en la operación escribimos algún número en forma decimal. corchetes y llaves Mathematica distingue entre paréntesis. ya que si confundimos su uso y escribimos por ejemplo Sin{x} o Sin(x) en lugar de Sin[x]. el programa nos lo hará saber mandándonos un mensaje de error de color azul./3 + 2/7 Out[5] = 0. Así.2 Paréntesis.619048. 4.Introducción a Mathematica lo obtendremos con un número 10 cifras signiﬁcativas. por ejemplo a la hora de calcular límites de funciones. al teclear In[5] := 1. • Los corchetes [·] se usan para escribir el argumento de una función bien sea matemática. Así el paréntesis de In[6] := (1 + 3)/7 4 Out[6] = 7 se usa para indicar que primero hacemos la suma 1+3 y luego dividimos entre 7. Mathematica distingue así entre operaciones algebraicas exactas y operaciones numéricas aproximadas. corchetes y llaves. y para escribir un número x real con seis cifras signiﬁcativas escribimos N[x]. En general es conveniente tener claro en qué momento se han de emplear los paréntesis. Asi por ejemplo. También se usan para construir conjuntos o listas de objetos matemáticos. Cada uno de estos elementos realiza una labor bien diferenciada en la estructura interna del programa. Por ejemplo la función sin x se escribe Sin[x]. si escribimos In[7] := 1 + 3/7 10 Out[7] = 7 vemos como el resultado cambia notablemente al realizarse en primer lugar la división y posteriormente la suma. los corchetes y las llaves. 65 .
22 (b) (2 − 3. x] = logb x Sin[x] = sin x Cos[x] = cos x ArcSin[x] = arcsin x ArcCos[x] = arccos x Tan[x] = tan x ArcTan[x] = arctan x n! = factorial de n Round[x] = parte entera de x Abs[x] = |x| FactorInteger[n] = f actores primos de n 66 .52 ¶56 1+i 3−i (f) 2 × 102 + 3 × 10−3 (i) (1 + i)7 (d) (3 + 4 ∗ 5) 5. si intentásemos calcular el determinante de una matriz no cuadrada.987 (e) µ 2.1 Errores Puede ocurrir que al teclear una operación en Mathematica y pulsar las teclas mayúscular + enter.1)23 (c) 3.Introducción a Mathematica Actividad 1 Realizar las siguientes operaciones con Mathematica: (a) 2. Esto ocurre cuando hay algún tipo de error o problema que el programa detecta. No debeis nunca de olvidar que el programa calcula rápidamente. Estos errores pueden ser básicamente de dos tipos: • Errores en la sintaxis de una sentencia.4+32 4∗7. Por ejemplo. el programa nos devuelva una salida conteniendo frases de color azul. el programa puede devolver un resultado erróneo aunque no nos escriba frases azules. Otras veces. 4. aunque esté bien escrita. 4. Por ésto es necesario saber qué estamos esperando de la operación que hemos pedido que el programa nos haga para así criticar el resultado y valorarlo en su justa medida. Es decir el programa no detecta ningún error a pesar de que éste existe.3 Funciones matemáticas de interés Mathematica posee una serie de funciones matemáticas predeﬁnidas que se escriben del siguiente modo: √ Sqrt[x] = x Exp[x] = ex Log[x] = log x Log[b.6 (g) (1 + i)(3 − i) (h) Actividad 2 Dar los resultados de la aplicación 1 con 15 cifras signiﬁcativas. • Errores producidos porque la expresión matemática o la operación realizada tiene algún problema. pero es tonto. Por ejemplo al escribir [1 + 2] ∗ 3 en vez de (1 + 2) ∗ 3 o N(π) en vez de N[π].75 + 8.2.3 ∗ 4 2 − 4.
respetando la sintaxis totalmente.4 Aprovechando cálculos anteriores A veces es posible que tengamos que hacer una sucesión de cálculos consecutivos de manera que cada nueva operación se base en la anterior.41421 N[Sqrt[7]. 67 . que nos proporciona un valor del número π con diez cifras decimales.34202. Actividad 3 Calcular los siguientes valores con 6 cifras signiﬁcativas: (a) sin 30◦ + cos 15◦ (b) log2 256 (c) e10! ¯ ¯ p 2 12 (f) ¯arcsin 1 + log 2 ¯ √ ¯ (d) (sin 1) (e) log 34 + e ¯ 2 2 4.6457513111. Parece necesaria entonces una sentencia que nos remita a resultados anteriores sin tener que escribir la operación de nuevo.Introducción a Mathematica Es importante destacar que hemos de escribir las funciones tal y como se detalla en la anterior tabla. Por ejemplo. para calcular el seno de 20 grados escribiríamos In[11] := Sin[20 Degree] //N Out[11] = 0. Si quisiéramos mayor precisión en las anteriores constantes. 10] 2.71828 ∞ √ i = −1 conversi´n de grados a radianes o Por ejemplo.10] Out[12] = 3.1415926536. Dicha sentencia es %. debemos escribir In[12] := N[Pi. Mathematica distingue entre letras mayúsculas y minúsculas. Entonces podemos calcular In[7] Out[7] In[8] Out[8] In[9] Out[9] In[10] Out[10] := = := = := = := = Sqrt[16] 4 Sqrt[2] √ 2 Sqrt[2] //N 1. y todas las funciones empiezan con letra mayúscula. Las constantes matemáticas más interesantes son: Pi E Inﬁnity I Degree = = = = = π ' 3.14159 e ' 2.
Esta operación podríamos hacerla del modo siguiente: In[13] Out[13] In[14] Out[14] := = := = Sin[20 Degree] //N 0. entonces 68 . Para remitirnos a un resultado obtenido dos pasos antes debemos escribir %%. Si x e y son dos variables que hemos deﬁnido con anterioridad. Entonces deberíamos hacer In[15] := G = 6. Algunas letras están asignadas ya por defecto por Mathematica y no pueden ser utilizadas para deﬁnir variables. hemos de tener en cuenta las siguientes reglas. Las variables pueden ser designadas por letras o por sucesiones de letras.34202 Cos[%] 0.5 Deﬁnición de variables y funciones Supongamos que tenemos que hacer una serie de cálculos en los cuales interviene repetidamente una constante. Supongamos que queremos deﬁnir la constante de la gravitación universal G = 6. 4. Una de ellas es N. el programa devolverá una sentencia azul de error. Aquí Cos[%] nos calcula el coseno del resultado obtenido en la salida 13. Para remitirnos a un resultado obtenido en la salida n podemos también escribir %n. y para resultados obtenidos k pasos antes escribiremos k símbolos %. la fuerza con la que se atraen dichos cuerpos se calcula In[16] := G ∗ 3ˆ2/(10ˆ2) Out[16] = 6. para después calcular el coseno de dicha cantidad.Introducción a Mathematica si queremos calcular cos (sin 20◦ ) tendríamos que calcular primero sin 20◦ . Es útil entonces deﬁnir variables con estos valores y funciones nuevas para minimizar el tiempo de escritura y hacer las operaciones de forma más ágil.003 × 10−12 . o al estudiar valores particulares de la función f (x) = (1. Si intentásemos escribir N=2.67 × 10−11 Si ahora tenemos dos cuerpos de masa 3 kilogramos separados a una distancia de 10 metros.67×10−11 con Mathematica. como por ejemplo la constante de la gravitación universal. Por ejemplo Cos[%13] también nos remite a la salida 13.67 × 10ˆ − 11 Out[15] = 6.942079.45)x + sin x. o bien In[18] := Clear[G]. Para desposeer a G de su valor debemos teclear In[17] := G = . A la hora de trabajar con variables en Mathematica.
• xy es una nueva variable formada por dos letras.. • xˆ2y es el producto x2 y y no x2y . y sólo ha escrito el último valor. y = 2. xn ] := expresión.45ˆx + Sin[x].45)x + sin x debemos escribir In[21] := f [x ] := 1. y = 2. 69 . Por ejemplo.”. obtenemos el valor de dicha función en 1. ha asignado el valor 1 a x. Entonces si tecleamos In[22] := Out[22] := f [1] 3. es decir la expresión In[20] := x = 1. x2 . para deﬁnir la función f (x) = (1. El símbolo _ se usa para indicar que la letra que lo antecede es una variable de la función. si en una misma línea queremos deﬁnir varias variables. z ] := {x ∗ y.. que como vemos es una función f : R3 → R2 . o escribir varias expresiones debemos separar estas con ”. no proporciona a continuación un Out[20] al pulsar mayúsculas + enter y asigna los mismos valores que en la sentencia anterior. Para eliminar la función debemos escribir In[23] := Clear[f.18975. z) = (xy. 2 a y y 3 a z. . Hemos visto como introducir funciones escalares. • 5x es el producto de 5 por x.Introducción a Mathematica • x y representará el producto x · y. x sin z). y . y. Por ejemplo In[19] := x = 1. z = x + y Out[19] = 3.. Por otra parte.. Para ello deﬁnimos la función como In[24] := f [x . Si al ﬁnal de la última expresión ponemos también ”. pero a menudo es necesario trabajar con funciones vectoriales.” la operación no proporciona ninguna salida. x ∗ Sin[z]}. por ejemplo la función f (x. Para deﬁnir nuevas funciones hemos de usar la siguiente estructura: nombrefunción[x1 . x] indicando tanto la variable como el nombre de la función. z = x + y.
especiﬁcando tanto la función como la variable respecto de la cual vamos a derivar. P i] {1. y) = sin(x + y) debemos escribir In[25] := D[Sin[x + y]. Para obtener las funciones coordenadas debemos escribir In[26] := Out[26] := o In[27] := Out[27] := f [1. Para calcular la derivada parcial con respecto a la variable y de la función f (x. x2 . P i][[1]] 1 obtenemos un vector del plano.. 70 . y] Out[25] = Cos[x + y]. x] ó D[f. 1. Así la segunda derivada de f (x) = sin x se calcula tecleando In[26] := D[Sin[x]. El comando que realiza ese cálculo con Mathematica es D[f.Introducción a Mathematica Si tecleamos ahora In[25] := Out[25] := f [1. 4. 0} f [1. {y.. Para calcular la derivada n—ésima de f (x) . {x. 1.6 Derivadas de funciones Supongamos que tenemos una función de una variable real f (x1 . n}]. 1.. xn ) a la que queremos calcular su derivada o derivada parcial respecto de alguna de sus variables. x] Out[24] = Cos[x]. hemos de proceder con el comando D[f. . P i][[2]] 0 según queramos trabajar con la primera o segunda coordenada. 3}] Out[27] = −Cos[x + y]. {x. xi ]. 2}] Out[26] = −Sin[x] y ∂3f ∂y3 de la función f (x. Por ejemplo si queremos calcular la derivada de f (x) = sin x escribiremos In[24] := D[Sin[x]. y) = sin(x + y) sería In[27] := D[Sin[x + y].
2 √ (c) y(x) = 1 + x2 . Así. . la variable de la función. 1 0. {x.Introducción a Mathematica Si ahora queremos calcular derivadas parciales de funciones respecto de diferentes variables hemos de indicarlo del modo siguiente D[f.5 -1 71 . donde indicamos la función. 2Pi}]. Así por ejemplo ∂ 2f de la función f (x. tenemos el comando Plot[f [x]. de la ecuación y 0 y = x. Exercise 1 Calcula las derivadas de las siguientes funciones: (a) f (x) = log (sin x) . x1 }].. de la ecuación y 00 + 2y 0 + y = ex . (d) y(x) = 1 (ex − 2xe−x − e−x ). Para ello hemos de darle tanto la función. (b) y(x) = 1 (e−x + 1). (b) f (x) = ex arcsin x . 2π] escribimos In[30] := Plot[Sin[x]. 0. y) = sin(x + y) se calcula escribiendo ∂x∂y In[28] := D[Sin[x + y]. xn ].5 1 2 3 4 5 6 -0. (a) y(x) = 2 − x + x2 . y un intervalo [x0 . para representar la función f (x) = sin x en el dominio [0. 4. (d) f (x) = xx . log4 (x2 +10) (c) f (x) = 1 + Exercise 2 Demostrar que las funciones siguientes satisfacen la ecuación diferencial que aparece a su lado.7 Representación gráﬁca de funciones Mathematica permite hacer representaciones gráﬁcas de funciones de una y varias variables.. x. x1 . de la ecuación y 0 + 2xy = x. y] Out[28] = −Sin[x + y]. como el dominio de deﬁnición de ésta. x0 . de la ecuación y 0 + y = x2 . x1 ] donde hacer la representación.. {x. x2 . Para la representación gráﬁca de funciones reales de variable real. 4 2 ³ x2 +tan 3x+ex √ x ´ .
{x.5 1 2 3 4 5 6 -0. fn [x]}. x1 }]. Para volver a representar gráﬁca una función ya representada previamente tenemos el comando Show[%n]. x0 . f2 [x]..5 1 2 3 4 5 6 -0. 2] . {x. 1 0. (d) f (x) = ex cos x en el dominio [−5. si escribimos In[33] := Show[%31]. Sin[2x]}.. en el dominio [−π. 0. ¡ 1+x ¢ (c) f (x) = sin 1−x2 en el dominio [−2. 2P i}] 1 0.Introducción a Mathematica Para representar varias funciones a la vez hemos de escribir todas las funciones que deseemos representar entre llaves y separadas por comas.. π] .5 -1 obtenemos una nueva representación gráﬁca simultánea de las funciones sin x y sin 2x. (b) f (x) = ex 1+x 1−x2 Exercise 4 Representar conjuntamente las gráﬁcas de los apartados (a). 5] . Así. 2] .5 -1 generaremos una representación gráﬁca simultánea de las funciones sin x y sin 2x. Si escribimos entonces In[31] := Plot[{Sin[x]. . Exercise 3 Representar gráﬁcamente las siguientes funciones de una variable: (a) f (x) = 1+x 1−x2 2 en el dominio [−2. (b) y (c) del ejercicio anterior. 72 . 2] . es decir Plot[{f1 [x]. (e) f (x) = ex cos x en el dominio [−2.
x]. x] Out[1] = {{y[x] → E 2 C[1]}} y la solución de la ecuación diferencial es de la forma y(x) = C[1]ex que proviene de la integración. indicando la ecuación y las variables dependiente e independiente. hemos de aprender a escribir ecuaciones diferenciales de manera que Mathematica las entienda.1 Ecuaciones diferenciales de primer orden Veamos cómo Mathematica es capaz de resolver ecuaciones diferenciales ordinarias. Resolución de ecuaciones diferenciales: orden uno y ecuaciones linales. 73 2 /2 x2 . En primer lugar. 5. donde C[1] es la constante . Así para resolver la ecuación y 0 = xy escribiremos In[1] := DSolve[y 0 [x] == x ∗ y[x].Capítulo 5 Ecuaciones diferenciales con Mathematica Sumario. y[x]. Esto se hace siguiendo la siguiente forma y 0 [x] == f [x. y[x]]. y) y(x0 ) = y0 . y) como problemas de condiciones iniciales de la forma ½ 0 y = f (x. Para calcular todas las soluciones de dicha ecuación diferencial tenemos la sentencia DSolve[y 0 [x] == f [x. y[x]. y[x]]. Podemos resolver tanto ecuaciones diferenciales de la forma y 0 = f (x.
calcular la población al cabo de 10 años. (e) y 0 = x − y. de agua pura.Ecuaciones diferenciales con Mathematica Para resolver problemas de condiciones tenemos que utilizar la sentencia anterior escribiendo la ecuación diferencial y la condición inicial entre llaves y separadas por comas. Exercise 5 Resolver las siguientes ecuaciones diferenciales de orden uno: (a) yy 0 = cos t. (f) 2xy 3 + 3x2 y 2 y 0 = 0. (d) Encontrar la concentración de sal en el tanque cuando éste tenga 25 l. la distancia entre (x. Si inicialmente había 10 personas que sabían la noticia y a los 3 dias la conocían 100 personas. y sale del tanque a razón de 2.6 l/min. de sal por litro entra en el tanque a razón de 1. Exercise 7 La población de medusas del Mar Menor varía de manera proporcional a la cantidad de medusas que hay en ese momento. Así el problema ½ 0 y = xy y(1) = 2 se resuelve escribiendo In[2] := DSolve[{y 0 [x] == x ∗ y[x]. y) de la misma. y) y el origen de coordenadas es igual a la longitud del segmento de la recta normal comprendido entre (x. Exercise 9 La velocidad a la que se transmite un noticia en un grupo es directamente proporcional al número de individuos que aun no la conocen. (c) y2 + 2yex + (y + ex )y 0 = 0. 2 (d) 2xy 3 + 3x2 y 2 y 0 = 0. x] 1 x2 1 x2 Out[2] = {{y[x] → 2E− 2 + 2 }} cuya solución es y(x) = 2e− 2 + 2 .000 individuos y al cabo de 2 años dicha población se triplicó. determinar cuanta gente lo sabrá al mes de producirse la noticia (tomar como población de España 40. y(0) = 0. Calcular la población de medusas para cada instante de tiempo t y calcula su límite cuando t → +∞. (b) Hallar la cantidad de agua en el tanque cuando la concentración de sal sea máxima. de agua. Si inicialmente la población de medusas era de 100. (c) Calcular la mayor cantidad de sal que llega a haber en el tanque en un momento dado. Se pide: (a) Determinar la concentración de sal en el tanque en cualquier tiempo.000). y[1] == 2}. y(π) = 3. Exercise 6 Hallar la familia de curvas que cumple que para todo punto (x. Una solución salina con 100 gr. y[x]. y(2) = 4. (b) y 0 = (1 + x)(1 + y). 74 .000. y) y el punto de corte de la recta normal con el eje x. En virtud del resultado obtenido ¿te parece acertado el modelo? ¿qué pegas le encuentras? Exercise 8 Un tanque contiene 40 l.3 l/min.
1) hemos de escribirla de la forma presentada en (5. Así. Por ejemplo.. (d) y 00 + 5y 0 + 6y = x3 cos x. (e) 2x2 y 00 + xy 0 + y = 0. de donde la solución general de la ecuación diferencial es y(x) = −x + C[1]e−x + C[2]ex . y[x]. En primer lugar. y 0 (0) = 1. Exercise 10 Resolver las siguientes ecuaciones diferenciales lineales con coeﬁcientes constantes: (a) 2y 00 + y 0 + y = x.Ecuaciones diferenciales con Mathematica 5. + pn−1 [x]y 0 [x] + pn [x]y[x] == f [x]. y 0 [0] == 0}. para resolver la ecuación y 00 − y = x In[1] := DSolve[y 00 [x] − y[x] == x. y[x]. para resolver el problema anterior ½ 00 y −y =x y(0) = 1... y(1) = 1. 1 ≤ i ≤ n. 0 [x] + . Finalizaremos estas prácticas indicando como resolver ecuaciones diferenciales lineales de orden mayor que uno.. donde C[1] y C[2] son dos constantes que provienen de la integración.2) hemos de escribir Nuevamente. Podemos además resolver problemas de condiciones iniciales escribiendo éstas entre llaves.. son funciones continuas. hemos de aprender a escribir la ecuación (5. y(0) = 0. + pn−1 (x)y 0 + pn (x)y = f (x). y 0 (0) = 0.2 Ecuaciones diferenciales lineales.. hemos de escribir lo siguiente: In[2] := DSolve[{y 00 [x] − y[x] == x. Por ejemplo.. x] obteniendo Notemos aquí que las condiciones iniciales C[1] y C[2] han desaparecido al aplicar las condiciones y(0) = 1 e y 0 (0) = 0. es decir.0 [x] + p1 [x]y 0 .2) dentro de una sentencia DSolve indicando la variable dependiente en primer lugar y la independiente a continuación.. x] que proporciona la siguiente respuesta Out[1] = {{y[x] → −x + E−x C[1] + Ex C[2]}}. Esto se hace escribiendo: y 0. (c) y 000 + 3y 00 + 3y 0 + y = ex . y 0 (1) = 0. (b) y 00 + 2y 0 + y = xex . y[0] == 1.1) donde f (x) y las funciones pi (x). para resolver la ecuación (5. n n−1 (5.1) de manera que el programa Mathematica la entienda. expresiones de la forma y n) + p1 (x)y n−1) + . 75 Out[2] = {{y[x] → (E−x x)(E2x − Ex x)}}. (5. . la sentencia que se utiliza para resolver ecuaciones de este tipo es DSolve.
3.3) es de la forma p p y(t) = c1 sin( k/mt) + c2 cos( k/mt). k = 0. k = 1 N/m. puede comprobarse que para cualquier condición inicial.1 Aplicaciones de las ecuaciones lineales con coeﬁcientes constantes. .4) p donde A recibe el nombre de amplitud. si suponemos el cuerpo libre de rozamiento y lo desplazamos verticalmente respecto de su posición de equilibrio. Se pide calcular las ecuaciones del movimiento para los siguientes valores de la masa y la constante de recuperación del muelle: (a) m = 1 kg. e y 0 (0) = 0.Ecuaciones diferenciales con Mathematica 5. Movimiento armónico simple.3) donde c1 y c2 son dos constantes reales que se calcularán una vez tengamos las condiciones iniciales y(0) e y 0 (0). (5. dicho cuerpo comienza a moverse según la ecuación diferencial my 00 + ky = 0. Exercise 11 Supongamos que desplazamos el cuerpo de la posición de equilibrio 1 m. e y 0 (0) = 1. 76 y(t) = A sin(ωt + ϕ). (5. (b) m = 2 kg. c2 = A sin ϕ.3 5. Si expresamos c1 y c2 en coordenadas polares ½ c1 = A cos ϕ. la solución de la ecuación (5. obtenemos la expresión donde m es la masa del objeto y k es la constante de recuperación del muelle. el movimiento armónico simple es aquél producido al colocar una masa en muelle como muestra la siguiente ﬁgura Entonces.5 N/m. ω = + k/m se conoce como frecuencia y ϕ como fase inicial. Como sabemos. Dado que la masa m y la constante k son positivas.
pedimos resolver el siguiente problema..5) my 00 + cy 0 + ky = 0. frecuencia y fase inicial de los movimientos anteriores. [0. En este caso se presentan tres situaciones perfectamente diferenciadas: • Movimiento sobreamortiguado. Dibujar las gráﬁcas de las funciones obtenidas al resolver las ecuaciones anteriores en el intervalo [0. Calcular las ecuaciones del movimiento en los siguientes casos: (a) m = 1 kg.. 5. c2 − 4mk > 0. (b) m = 1 kg. c2 − 4mk = 0.3. Exercise 13 Hacer un estudio comparativo de los movimientos generados por los apartados de la Actividad 12 con el de un movimiento armónico simple tal que m y k tengan el valor que tenían en el apartado y c = 0. 10π] y comprobar que son periódicas. ¿Qué información puedes obtener de las gráﬁcas de las funciones? Caracteriza cualitativamente el movimiento sobreamortiguado. Como orientación podeis tomar los intervalos dados en el ejercicio 12 para hacer las comparaciones. Exercise 12 Desplazamos de nuevo el cuerpo de la posición de equilibrio 1 m. 20π]. • Movimiento subamortiguado. Si suponemos que sobre el cuerpo colgado del muelle anterior se ejerce una fuerza de rozamiento proporcional a la velocidad con constante de proporcionalidad c. 2π]. críticamente amortiguado y subamortiguado. • Movimiento críticamente amortiguado. c = 3 N · sg/m y k = 1 N/m.. (c) m = 1 kg. Determinar qué tipo de movimiento se obtiene en cada caso y hacer sucesivas representaciones gráﬁcas de las funciones anteriores en los intervalos [0. e y 0 (0) = 2. (d) m = 1 kg. 8π] y [0. sabemos que la ecuación del movimiento se escribe de la forma (5. (e) m = 1 kg. 4π]. 77 . c = 4 N · sg/m y k = 4 N/m. Representar a la vez las gráﬁcas de las funciones obtenidas al resolver la ecuación amortiguada y la no amortiguada. calculando el periodo de éstas. Obtener además la amplitud. c = 2 N · sg/m y k = 1 N/m. [0. c = 3 N · sg/m y k = 5 N/m..2 Movimiento amortiguado. c = 1 N · sg/m y k = 4 N/m. En esta sección. c2 − 4mk < 0.Ecuaciones diferenciales con Mathematica (c) m = 1 kg. k = 4 N/m..
(b) F (t) = cos t. (c) F (t) = sin t. (b) y 00 + 2y 0 + 2y = e−t . [0. c = 2 y k=2 y 00 + 2y 0 + 2y = 0. (d) F (t) = cos(2t). (h) F (t) = sin(2. una resistencia R.4 Aplicación a los circuitos eléctricos Consideremos un circuito eléctrico que lleve en serie una bobina de inductancia L. Supondremos ahora que el cuerpo que inicialmente considerábamos colgado de un muelle. 10].e. (e) F (t) = sin(2t). (g) F (t) = t.m. V (t).3.3 Movimiento forzado. (i) F (t) = cos(1. (f) F (t) = et .1t). Hablaremos entonces de movimiento forzado. un condensador de capacidad C y que es alimentado por una f. [0. x(0) = 1. además de estar sujeto a fuerzas de rozamiento. está afectado por una fuerza F (t) que modiﬁca y condiciona su movimiento. Obtener a continuación para cada uno de los siguientes apartados la solución particular de la ecuación no homogénea correspondiente y la representación gráﬁca por separado y conjunta de la solución particular de la ecuación no homogénea y de la solución de los problemas de condiciones iniciales correspondientes en los intervalos [0. cuya ecuación de movimiento viene dada por la ecuación diferencial my 00 + cy 0 + ky = F (t). Hacer un estudio de las gráﬁcas de las funciones resultantes en un intervalo de la forma [0. ¿Qué gráﬁcas parecen acotadas? ¿Qué explicación le das para el caso en que las gráﬁcas no estén acotadas? Exercise 15 Resolver la ecuación diferencial correspondiente a un muelle donde m = 1. 5] y [0. (c) y 00 + 2y 0 + 2y = t2 + 1. x(0) = 1. x0 (0) = 0. aumentando el valor de T progresivamente. 2]. como en el ejercicio 12. según muestra la siguiente 78 . ¿Qué consecuencias puedes sacar a partir de las gráﬁcas obtenidas? 5.6) Exercise 14 Resolver el problema de condiciones iniciales ½ 00 y + 4y = F (t) y(0) = y 0 (0) = 0 en los siguientes casos: (a) F (t) = 1. 1].9t).Ecuaciones diferenciales con Mathematica 5. x0 (0) = 0. x(0) = 1. (a) y 00 + 2y 0 + 2y = cos t. (5. T ]. x0 (0) = 0.
R y C son constantes.8) Como puede apreciarse.7) y (5. V (t) = VC + VR + VL donde VC .Ecuaciones diferenciales con Mathematica ﬁgura Suponiendo que L. obtenemos la ecuación lineal de orden dos Lq 00 (t) + Rq0 (t) + q(t)/C = V (t). mediante física elemental se sabe que el voltaje generado V (t) se consume en todos los elementos del circuito. cabe el mismo análisis para circuitos que hicimos en el apartado anterior. Sabiendo que q(t) VC = . es decir.7) Teniendo en cuenta que la intensidad i(t) se deﬁne como la derivada de la carga q(t) obtenemos la ecuación en términos de la intensidad Li00 (t) + Ri0 (t) + i(t)/C = V 0 (t) (5. Exercise 16 Consideremos el circuito eléctrico de la ﬁgura. C donde q(t) es la carga en cada instante de tiempo. las ecuaciones (5.8) son idénticas a la ecuación que proviene de la vibración de un muelle. VR = Rq 0 (t) y VL = Lq00 (t). (5. 79 . la resitencia y la bobina respectivamente. Así. VR y VL representan la diferencia de potencial entre el condensador.
L = 2H. R = 2Ω. (b) C = 1F . V (t) = sin t. R = 3Ω. L = 0H. 80 .25F . V (t) = −t cos t. suponiendo que el circuito está descargado (i(0) = i0 (0) = 0): (a) C = 1F . V (t) = e3t . R = 4Ω. V (t) = sin t. L = 1H. L = 1H. V (t) = t2 . haciendo un estudio gráﬁco de la misma.Ecuaciones diferenciales con Mathematica Calcular la intensidad de corriente que pasa por los cables de dicho circuito en los siguientes casos. V (t) = et cos(2t). (c) C = 2F . L = 1H. R = 1Ω. L = 0H. (d) C = 1F .5F . R = 1Ω. (f) C = 0. (e) C = 0.
A continuación vamos a explicar algunas sentencias que se utilizan para hacer programas con Mathematica.5. si queremos asignar a la letra q el valor 6.1. escribiremos en la pantalla q = 6.1 6.1 Operaciones y deﬁnición de variables. vamos a contestar a la primera pregunta: para interrumpir un programa o una evaluación de Mathematica basta con pulsar las teclas Alt + coma. En primer lugar. se necesita una explicación más amplia que a continuación desarrollamos. Las principales herramientas que vamos a necesitar para hacer los programas son las siguientes: • ¿Cómo establecer comunicación entre el ordenador y quien ejecuta el programa? • ¿Cómo contruir bucles con Mathematica? • ¿Cómo deﬁnir funciones y variables? • ¿Cómo interrumpir la ejecución de un programa cuando éste esté bloqueado? Vamos a explicar como hacer cada una de estas cosas.Capítulo 6 Programación en Mathematica Sumario. Por ejemplo. Para las otras cuestiones. 6.5 los que producirá la siguiente respuesta en la pantalla del ordenador In[1] := q = 6. para deﬁnir una variable numérica basta con tomar una letra o conjunto de letras e igualarlas al valor deseado. Representación gráﬁca de datos discretos.5 Out[1] = 6.5 81 . para poder construir nuestro programa. Generalidades de la programación en Mathematica. Deﬁnición de variables y funciones Como sabemos.
Programación en Mathematica Si lo que queremos es deﬁnir una función. Podemos deﬁnir una variable vectorial. para deﬁnir la función f (x) = ex cos x hemos de escribir In[2] := f[x ] = Exp[x] ∗ Cos[x] Out[2] = Exp[x] ∗ Cos[x]. Escribiendo Clear[q] desposeemos a q del valor 6. 9} dejando el valor de a[2] vacio. Para borrar las variables o las funciones utilizaremos la sentencia Clear[var o func]. o Así. hemos de hacerlo indicando la variable entre corchetes según el siguiente esquema: Nombre[var ] = funci´n. 9}. por ejemplo de dimensión tres haciendo la siguiente operación In[3] := a[1] = 5 Out[3] = 5 In[4] := a[2] = 2 Out[4] = 2 In[5] := a[1] = 9 Out[5] = 9 y habremos introducido el vector (5. 9). a[2]. 3}] Out[7] = {5. 2. Podemos dejar sin asignar valores al vector.5 asignado anteriormente. Otro tipo de variables bastante interesantes son las vectoriales. Si queremos saber el valor de la variable a hemos de escribir In[6] := ?a Global 0 a a[1] = 5 a[2] = 2 a[3] = 9 o bien In[7] := Table[a[i]. 3}] Out[7] = {5. {i. 82 . Si hacemos In[8] := Clear[a[2]] y In[7] := Table[a[i]. {i. 2.
Si escribimos 7 > 4 y 26= 3 aparecerá In[2] := 7 > 4 && 2 ! = 3 Out[2] = True indicando que las dos sentencias son verdaderas. Podemos combinar esta operaciones para obtener otros tipos de operaciones y estructuras lógicas muy usadas en programación. Si ponemos una verdadera y una falsa y escribimos In[3] := 7 > 4 && 2 = 3 Out[3] = False mientras que In[4] := 7 > 4 || 2 = 3 Out[4] = True 83 . aparecerá en pantalla In[1] := 10 < 8 Out[1] = False indicando la falsedad de la aﬁrmación.2 Operaciones + − ∗ / Suma Resta Multiplicación División Como sabemos las cuatro operaciones básicas que tenemos son junto con el símbolo =. si escribimos 10 < 8. Como las más interesantes resumimos las siguientes: Operaciones de Test == Igual ! = Distinto < Menor que > Mayor que <= Menor o igual que >= Mayor o igual que Operaciones Lógicas ! o Not Negación && Y || O True Verdad False Falsedad Por ejemplo.Programación en Mathematica 6.1.
Si escribimos In[6] := x+ = 2 Out[6] = 9 el valor de x será 9. Otras combinaciones de símbolos que pueden resultar útiles son las siguientes: x++ x−− ++x −−x x+ = k x− = k x∗ = k x/ = k Aumenta el valor de x una unidad Disminuye el valor de x una unidad Preincrementa x una unidad Predisminuye x una unidad Aumenta el valor de x k unidades Disminuye el valor de x k unidades Multiplica x por k Divide x por k Por ejemplo. Basic o C. si hacemos x = 7 y realizamos las siguientes operaciones.2. los cuales han de ejecutarse un número de terminado de veces. explicamos a continuación una serie de sentencias útiles. Recuerdan en gran medida a expresiones de lenguajes de programación como el Pascal. 6. De un modo parecido funciona con −−.1 If La sentencia ”If” se emplea para establecer elementos condicionantes en un programa. If[condici´n. que es verdadero. obtenemos In[6] := x + + Out[6] = 7 donde el valor de x ahora es 8. igual que el que aparece en pantalla. o 84 .2 Construcción de programas con Mathematica. f ]. pero en pantalla se escribe su antiguo valor. Las diferentes sintaxis con las que puede manejarse son las siguientes. Si hubiéramos escrito In[6] := + + x Out[6] = 8 ahora el valor de x sigue siendo 8 y el programa escribe el valor actual. Uno de los temas centrales en los programas que vamos a realizar es la construcción de bucles o procesos iterativos. 6.Programación en Mathematica Finalmente In[5] := Not[7 < 4] Out[5] = True dado que lo contrario de 7 < 4 es 7 > 4. Para la una correcta creación de bucles.
b. 3. 2 + 1. y si es falsa no hará nada. 2 + 2] Out[1] = 4 pero si escribimos In[2] := If[2 > 1. x > 5. 6. 0]. 2 + 2] no se obtendrá ninguna salida. Evalúa la condición y la compara con C1 . 1. f2 ]. 6. Es otra sentencia muy utilizada para deﬁnir funciones a trozos. si la condición es cierta el programa ejecutará f1 y si es falsa ejecutará f2 .3 Switch Esta sentencia también es útil para producir toma de decisiones en un programa. Es útil a la hora de deﬁnir funciones a trozos. 2 + 2] Out[3] = 4.. . c] Out[1] = c. compara la condición con C2 y si ambas coinciden ejecuta f2 y así sucesivamente. o Ahora. 0. C2 . f2 . C1 .. In[1] := If[2 > 1. Su estructura es la siguiente Switch[condici´n. deﬁnirá la función ⎧ 2 ⎨ x si x < 0 x3 si x > 5 h(x) = ⎩ 0 si 0 ≤ x ≤ 5. f1 . f2 . xˆ2.2 Which La sentencia ”which” también se emplea para introducir elementos de toma de decisiones en un programa. . Por ejemplo In[1] := Switch[3. f1 .]. test2 . Si ambas coinciden ejecuta f1 . A continuación. Otro tipo de expresión es la siguiente If[condici´n.Programación en Mathematica Si la condición es verdadera el programa ejecutará f . Por ejemplo.2. o y funciona de la siguiente manera. a... Por ejemplo In[1] := h[x ] = Which[x < 0. f1 . La estructura de la sentencia es la siguiente Which[test1 .2. xˆ3. True. Por ejemplo In[3] := If[1 > 2.]. 85 .
Para obtener la suma de los cubos de los números pares entre el uno y el 100 debemos escribir In[2] := s = 0. por lo que hemos de indicarle al ﬁnal que lo escriba. la sentencia While no escribe nada en la pantalla. Nótese que. o Do[expresi´n. (6. En la segunda sentencia (6.3) evaluamos la expresión desde un valor mínimo que nosotros escogemos.5 While Otra manera de construir bucles es empleando la sentencia While. While[(i = i + 1) <= 10. pero ahora aumentamos el valor de i en cada paso de di en di unidades. s Out[1] = 385.1) hemos de escribir In[1] := s = 0.4) La sentencia (6.2. 10. Do[s = s + iˆ3. 0. Para ello necesitamos en primer lugar saber como se comporta y cual es su sintaxis. La sintaxis es la siguiente: While[test. i = 0. Es importante también darse cuenta de como escribir la condición y el papel que juegan los paréntesis en la misma. aumentando su valor de uno en uno. como por ejemplo calcular la suma 10 X i2 . {i. 6. 2. Es por ello que al ﬁnal escribimos s en In[1].Programación en Mathematica 6. di}].4 Do Esta sentencia es útil para la construcción de bucles o procesos repetitivos. Son las siguientes: Do[expresi´n.2) evalúa la expresión desde el valor i = 1 hasta el valor i = imax .1) i=0 Veamos como hacerla empleando esta sentencia. Do[s = s + iˆ2. expresi´n]. imin . 86 . para obtener la suma de (6.2. Por último. al igual que pasaba con Do. o Do[expresi´n.2) (6.1) podríamos escribir In[1] := s = 0. A diferencia de la anterior. Hemos de remarcar que la sentencia Do no escribe nada en pantalla. 100. s Out[2] = 13005000. 2}]. Así. imax }]. {i. imax . hasta un valor máximo i = imax .3) (6. Por ejemplo. o que produce la ejecución de la expresión mientras es test de un resultado positivo. para obtener el valor de la suma (6.4) ejecuta la expresión i veces. esta sentencia incorpora la posibilidad de toma de decisiones. o (6. s Out[1] = 385. {i. {i}]. {i. 1}]. (6. i = imin . sólo ejecuta la operación que se le indica. s = s + iˆ2].
Así. i + +. i < 0. If[s > 400. i + +. s Out[1] = 385. s Out[1] = 385 que nos vuelve a calcular la suma (6. 6.1) hemos de escribir In[1] := s = 0. si escribimos For[i = 1. expresi´n]. i < 0. i = 0. i + +. si escribimos For[i = 1. si escribimos In[1] := s = 0. Por ejemplo. Break[ ]]. Por ejemplo. que producirá una salida del bucle al superar s el valor 400. 87 . For[i = 1. Goto[hola]]. i <= 10. Esto provocaría una situación en la cual el ordenador se ”cuelga”. Nuevamente no escribe nada en pantalla. Label Estas dos sentencias se emplean para conectar dos sitios determinados de un programa. para intentar escribir bucles que no seamos capaces de construir usando las sentencias anteriores. Estas sentencias han de emplearse en último lugar. Break[ ]]. s = s + iˆ2]. por lo que para obtener la suma (6. Escribiendolo en pantalla tendríamos In[1] := s = 0. If[s > 400.1).6 For Otro modo de construir un bucle combinado con toma de decisiones es la sentencia For. La sentencia Break puede ser empleada para impedir este proceso.2. s = s + iˆ2. i + +. i + +.8 Goto. La sintaxis de ambas es Label[etiqueta]. incremento.Programación en Mathematica 6. For[i = 0. o que partiendo de un valor inicial. s = s + iˆ2]. i < 0. Goto[etiqueta]. s = s + iˆ2].2. teniendo la sintaxis Break[ ].2. 6. s = s + iˆ2] sumaríamos indeﬁnidamente ya que i nunca va a ser menor que cero. evalúa la expresión a la vez que incrementa el valor inicial hasta que el test falla. test. s Out[1] = 0. Label[hola]. If[i <= 10.7 Break Frecuentemente al trabajar con bucles puede llegarse a una situación que provoca la repetición indeﬁnida del mismo. cuya sintaxis es la siguiente: For[inicio.
Por ejemplo.10 Input La sentencia Input se emplea para introducir datos al programa desde el exterior del mismo cuando éste está en ejecución. Se utiliza para establecer un diálogo entre el programa y el usuario. hay sentencias que no escriben nada en la pantalla. Como hemos visto. Cuando se escribe (6. utilizando la expresión Print[”F rase”]. Con Print podemos conocer el valor de las evaluaciones producidas. i = 0. La sintaxis es la siguiente: Input[ ] o Input[Anotaci´n].5) aparecerá un cuadro de diálogo que habrá que rellenar y que suelen ser datos útiles del programa. While[(i = i + 1) <= 4.2. Puede suceder que al escribir frases con la sentencia Print el programa cambie las palabras de orden. pero con la anotación que hayamos decidido poner. podemos ir viendo paso a paso cómo se calcula el valor de la suma 4 X i=0 i2 escribiendo In[1] := s = 0. Print[s]] 1 5 14 30. por ejemplo.9 Print Esta sentencia se utiliza exclusivamente para escribir en la pantalla. Así. se escribimos Input[escribe tu edad] 88 .6) (6. Para evitar esto basta con escribir la frase entre comillas. o (6.6) aparecerá el mismo cuacdro. s = s + iˆ2.5) Cuando se escribe (6.2.Programación en Mathematica 6. 6.
y1 = y0 + h ∗ f[t0. t1 = t0 + h.3 Un programa Vamos a diseñar un programa que aproxime numéricamente mediante el método de Euler la solución del problema de condiciones iniciales ½ 0 y = ty. y0 = Input[“Valor de inicial de la funci´n”]. n = Input[“N´mero de pasos”]. por ejemplo 22. aparecerá en pantalla lo siguiente: In[1] := Input[escribe tu edad] Out[1] = 22. y ] := t ∗ y. t0 = t1]. f[t . 89 . i < n + 1. donde las condiciones iniciales t0 e y0 . o tf = Input[“Tiempo ﬁnal”]. For[i = 1. t0 = Input[“Tiempo inicial”]. el instante ﬁnal tf en el queremos conocer el valor de la solución. i + +. y0 = y1. u h = (tf − t0)/n. y aceptado. Print[“El valor aproximado es ” y0].Programación en Mathematica aparecerá un cuadro de la siguiente forma que habrá que cumplimentar y aceptar para que desaparezca. y el número de pasos con el que este valor se aproximada han de introducirse tras la ejecución del programa. y0]. y(t0 ) = y0 . 6. Una vez introducido el valor.
Programación en Mathematica
Presentaciones gráﬁcas
A menudo, aparte de los datos numéricos en sí, es interesante tener una representación gráﬁca de los mismos. En el caso de la aproximación numérica de soluciones de ecuaciones diferenciales, éstos son dados en forma de sucesión de números reales yn , con n variando de 1 a un valor N, de manera que yN es la aproximación a la solución para el tiempo tf que hemos elegido. En el caso de sistemas de ecuaciones diferenciales, tendremos una sucesión de vectores de los cuales, cada coordenada se corresponde con la aproximación de la función solución en el instante de tiempo dado. En cualquiera de los dos casos, es conveniente tener una manera de presentar los datos obtenidos mediante la integración numérica de una manera gráﬁca. En general, el comando ListPlot permite hacer tales presentaciones mediante la estructura, pero previamente, tenemos que tener una manera de introducir en Mathematica las sucesiones con las que deseamos trabajar. Dada una sucesión de números reales yn , ésta se puede introducir en Mathematica de dos maneras. La primera es hacerlo como una lista en una variable, que se haría con la sintaxis In[1] := Y = {y1 , y2 , ..., yn } Out[1] = {y1 , y2 , ..., yn } No obstante, si tenemos muchos datos, digamos 1000, no es operativo introducirlos de esta manera. Para ello tenemos la sentencia Append, que tiene la sintaxis Append[Y, yn+1 ] que agrega el dato yn+1 al ﬁnal de la lista Y . Veamos cómo: In[1] Out[1] In[2] Out[2] := = := = Y = {2, 5, 6} {2, 5, 6} Append[Y, 8] {2, 5, 6, 8}
Supongamos ahora que deseamos introducir los mil primeros términos de la sucesión dada por la recurrencia yn+1 = yn + 2, donde y1 = 1, es decir, yn = (1, 3, 5, 7, ...). Escribimos el siguiente programa Y = {1}; x0 = 1; For[i = 1, i < 1000, i + +, x1 = x0 + 2; Z = Append[Y, x1]; x0 = x1; Y = Z] Una vez ejecutamos el programa, almacenaremos en Y el valor de la sucesión. Una manera alternativa de introducir dicha sucesión es mediante los arrays, que en dimensión uno se escriben como Array[Y, n]. Esta sentencia genera automáticamente un vector que tiene por coordenadas Y [1],Y [2],...,Y [n], que en principio estará vacío de contenido. Para introducir la sucesión anterior escribiríamos el programa Array[Y, 1000]; x0 = 1; For[i = 1, i < 1001, i + +, Y [i] = x0; x1 = x0 + 2; x0 = x1] 90
Figura~6.1: Representación gráﬁca proporcionada por la salida del comando ListPlot. Si deseamos que en el array anterior el índice empiece a contar de cero, debemos usar la forma Array[Y, n, 0]. Diferentes alternativas para estas sentencias pueden consultarse en la ayuda del programa. Pasemos a ver cómo representar gráﬁcamente la información con el comando ListPlot, que tiene cualquiera de las siguientes estructuras ListPlot[Y ] ó ListPlot[Array[Y, 1000]] según hallamos introducido la sucesión de una u otra manera. Así, tanto el programa Y = {1}; x0 = 1; For[i = 1, i < 1000, i + +, x1 = x0 + 2; Z = Append[Y, x1]; x0 = x1; Y = Z]; ListPlot[Y ] como su alternativa Array[Y, 1000]; x0 = 1; For[i = 1, i < 1001, i + +, Y [i] = x0; x1 = x0 + 2; x0 = x1]; ListPlot[Array[Y, 1000]] nos devuelven la gráﬁca dada en la ﬁgura 6.1.
El comando ListPlot tiene la opción de unir dos puntos consecutivos con una recta. Tendría la sintaxis ListPlot[Y, PlotJoined− > True] ó ListPlot[Array[Y, n], PlotJoined− > True] Veamos como obtener una información gráﬁca de de aproximación que obtuvimos del problema de condiciones iniciales generado a partir de la ecuación y 0 = ty. Para ello, basta modiﬁcar el programa 91
Figura~6.2: Aproximación de la solución del problema de condiciones inciales. anterior según se indica f[t , y ] := t ∗ y; t0 = Input[“Tiempo inicial”]; y0 = Input[“Valor de inicial de la funci´n”]; o tf = Input[“Tiempo ﬁnal”]; n = Input[“N´mero de pasos”]; Array[s, n + 1, 0] u h = (tf − t0)/n; s[0] = y0; For[i = 1, i < n + 1, i + +, s[i] = y0 + h ∗ f[t0, s[i − 1]]; t1 = t0 + h; t0 = t1]; ListPlot[Array[s, n, 0], PlotJoined− > True]; obtenemos la ﬁgura 6.2 como aproximación de la solución con condición inicial y(0) = 1, donde el tiempo ﬁnal es 10 y el número de pasos 100. Como puede apreciarse en la gráﬁca 6.2, en el eje x aparecen datos hasta 100, que se corresponden con el número de pasos. Si hubiésemos elegido 1000 pasos, aparecería obviamente 1000 al ser éstos en número de elementos en la sucesión generada. Podemos hacer una representación gráﬁca donde el eje x recorra el intervalo de deﬁnición de la función solución de la siguiente manera f[t , y ] := t ∗ y; t0 = Input[“Tiempo inicial”]; y0 = Input[“Valor de inicial de la funci´n”]; o tf = Input[“Tiempo ﬁnal”]; n = Input[“N´mero de pasos”]; Array[s, n + 1, 0] u h = (tf − t0)/n; s[0] = y0; r = {{t0, s[0]}} For[i = 1, i < n + 1, i + +, s[i] = y0 + h ∗ f[t0, s[i − 1]]; t1 = t0 + h; t0 = t1 ; r1 = Append[r, {t0, s[i]}]; r = r1]; ListPlot[r, PlotJoined− > True]; es decir, deﬁniendo una sucesión r de vectores del plano, y representando ésta como muestra la ﬁgura 6.3 92
es necesario tener que representar a la vez dos sucesiones. g[t0]}]. Por ejemplo. Supongamos que queremos ver una representación gráﬁca conjunta de dicha función y su aproximación en el intervalo [0. para representar ambas funciones a la vez hemos de activar un paquete especial de Mathematica tecleando << Graphics‘MultipleListPlot‘ A continuación utilizamos la sentencia MultipleListPlot.3: Representación de la solución en el dominio de deﬁnición.. i + +.4. h = (tf − t0)/n. Sin embargo. las soluciones aproximadas y exacta de una ecuación diferencial para conocer o estimar la bondad de la aproximación. Otras opciones para esta sentencia pueden verse en el manual del programa. graf1 = Append[graf. la solución del problema ½ 0 y = ty. n = 100. y(0) = 1. generamos los valores de la solución exacta en los 100 puntos donde obtenemos la aproximación numérica con el programa g[t ] := Exp[tˆ2/2]. t1 = t0 + h. t0 = t1. t0 = 0. graf = {{t0. 1] mediante 100 pasos. graf = graf 1] Igualmente. F or[i = 1. graf}. {t0. A veces.Programación en Mathematica Figura~6. 93 . Por ejemplo. tf = 1. que permite representar varias sucesiones con la siguiente sintaxis MultipleListPlot[{r. i < 101. g[t0]}}. es y(t) = et 2 /2 . PlotJoined− > True] que da lugar a la gráﬁca 6. Para ello. ejecutamos el programa anterior para obtener la solución aproximada almacenada en la variable r.
Programación en Mathematica Figura~6. 94 .4: Representación conjunta de las soluciones aproximada y exacta del problema de condiciones iniciales.
existe una gran diferencia entre las ecuaciones anteriores. consideraremos ecuaciones en diferencias de la forma yn+k = f (n. son de órdenes 2 y 3. Por ejemplo. xn .Apéndice A Ecuaciones en diferencias Sumario. Solución de una ecuación en diferencias... Estabilidad de ecuaciones en diferencias lineales. nyn+3 − eyn+3 yn = yn+1 . Ecuaciones lineales: teoría general.. que llamaremos resueltas respecto de el mayor término de la sucesión yn . yn . no se va a poder despejar explícitamente el término yn+3 . .. A. yn+k−1 ).. yn+k ) = 0. yn+1 . quedando la ecuación yn+2 = yn .1) entenderemos una solución xn de números reales de manera que veriﬁque xn+k = f (n. El número k recibe el nombre de orden de la ecuación. Aparte del orden. Ecuaciones en diferencias..1 Introducción F (n. es decir. Estabilidad local de ecuaciones en diferencias no lineales. Nosotros vamos a centrarnos en el primer tipo de ecuaciones.. xn+1 . las ecuaciones yn+2 − yn = 0. yn+1 . Transformada Z.1) siendo f : Λ ⊆ Rk → R una función. Por una solución de la ecuación (A. . A partir de este momento.. 95 . mientras que en la segunda ecuación tal operación no puede realizarse. (A.. . En la primera se puede despejar el término yn+2 . respectivamente. yn . xn+k−1 ). Una ecuación en diferencias es una expresión de la forma donde F : Ω ⊆ Rk+2 → R es una función deﬁnida sobre un subconjunto Ω de Rk+1 .
retarda la señal o sucesión de entrada una 96 . .. Por ejemplo El segundo elemento. y2 = 1. por ejemplo la sucesión constante xn = 1 es solución de la ecuación yn+2 = yn . Como vemos.. y1 = −1. También lo es la sucesión xn = (−1)n . En el caso de que las sucesiones a1 . Así.Ecuaciones en diferencias Así. ak sean n n n n constantes. n n donde a1 . que son de la forma yn+k + a1 yn+k−1 + . + ak yn = bn . son ecuaciones en diferencia lineales. siendo la primera no homogénea y la segunda homogénea y con coeﬁcientes constantes. ai = ai para todo n ≥ 0 y para todo i ∈ {1. una D dentro de un rectángulo. yn+2 − yn+1 − yn = 0. Llamaremos a estos problemas de condiciones iniciales. . y no homogéneas en caso contrario. bn son sucesiones de números reales. Este dispositivo genera una sucesión de salida yk para una sucesión de entrada xk de la siguiente manera. pero ésta es única si imponemos una serie de k condiciones iniciales. ak . El elemento marcado con una a dentro de un círculo ampliﬁca el dato de entrada la magnitud a ∈ R.... Dentro de las ecuaciones en diferencias. Veámoslo con el siguiente ejemplo que proviene de la electrónica (ver [Oga2]). consideremos el siguiente ejemplo.. una ecuación puede tener distintas soluciones.. la ecuación lineal se dirá de n coeﬁcientes constantes. . En general. Las ecuaciones yn+3 + nyn+1 − yn = 1.. Las ecuaciones lineales juegan un importante papel en la modelización de circuitos digitales... por su analogía con las ecuaciones diferenciales ordinarias. k}. Para ﬁjar ideas. tienen un espcial interés las llamadas ecuaciones lineales. esto es.. también distinguiremos entre ecuaciones homogéneas si bn = 0 para todo n ≥ 0. xn = (−1)n es la única solución de la ecuación ½ yn+2 = yn .
En el proceso. Combinando (A. (A. el elemento marcado con un símbolo S dentro de un círculo. se obtiene una ecuación de orden dos. construimos los llamados circuitos digitales. los datos intermedios rk vienen dados por la expresión rk = xk − ayk . Si complicamos el dispositivo.2) y (A. Aquí yk+1 = vk . Así Finalmente.Ecuaciones en diferencias unidad temporal.2) donde a es un número real.3) obtenemos la ecuación en diferencias de orden uno yk+1 + ayk = xk . Los datos de entrada vienen dados por la sucesión xk y los de salida por yk+1 = rk . 97 . como se muestra en la ﬁgura. Ésta representa uno de los tipos más sencillos de retroalimentación de una señal.3) (A. suma los datos que le llegan: Combinando varios de estos elementos. como el de la ﬁgura anterior.
). .4) zn n=0 P Nótese que (A. Propiedades básicas. Veremos a continuación cómo abordar el estudio de estas ecuaciones. 0. • Linealidad. Dadas las sucesiones xk e yk y α. y1 = 1.2 A. 0. n=1 Por ejemplo. +∞) = {z ∈ C : |z| > r} 1 z = z .) entonces su transformada Z es Z[δ](z) = 1 deﬁnida en todo el plano complejo. 1. β ∈ C. 0. z−1 . r. .1 Transformada Z Deﬁnición y propiedades básicas Consideremos una sucesión de números complejos xk . y0 = 0. 1. 98 ∞ X 1 1 = n z 1− n=0 A(0. de donde se obtiene la ecuación yk+2 + ayk+1 − byk = xk .. Si xk = (1.. A..4) es una serie de Laurent con parte regular x0 y parte singular ∞ xn z −n . Por ejemplo supongamos la ecuación ½ yk+2 + yk+1 − 2yk = 0. si δ = (1. y que n=1 por tanto convergerá en un disco de convergencia de la forma P donde r es el radio de convergencia de la serie de potencias ∞ xn z n . Se deﬁne la transformada Z de la misma como la serie ∞ X xn Z[xk ](z) = .2.. se veriﬁca Z[αxk + βyk ](z) = αZ[xk ](z) + βZ[yk ](z) para todo z en el dominio de deﬁnición de Z[xk ](z) y Z[yk ](z).Ecuaciones en diferencias vk+1 = rk . entonces Z[1](z) = siempre que |z| > 1. rk = xk + byk − avk . (A.
zn • Dada la sucesión xk y a ∈ C \ {0}. se tiene que Z[2 ](z) = k ∞ X 2n n=0 zn = 1 1− 2 z = z . Basta calcular Z[αxk + βyk ](z) = ∞ X αxn + βyn n=0 ∞ X n=0 zn = α ∞ X yn xn +β = αZ[xk ](z) + βZ[yk ](z). si xk = (1. tenemos la fórmula Z[xk+k0 ](z) = z k0 Z[xk ](z) − Demostración. deﬁnimos la nueva sucesión yk = xk+1 . Entonces Z[yk ](z) = Z[xk+1 ](z) = zZ[xk ](z) − zx0 . z−2 99 . 23 .).. zn = z = z ∞ X xn xn+1 =z z n+1 zn n=1 xn − zx0 = zZ[xk ](z) − zx0 . Dmostración. zn zn n=0 • Dada la sucesión xk .Ecuaciones en diferencias Demostración. (z/a)n Por ejemplo. 2. si k0 ∈ N y deﬁnimos yk = xk+k0 . Calculamos Z[a xk ](z) = k ∞ X an xn n=0 zn = ∞ X n=0 xn = Z[xk ](z/a). . En general. se veriﬁca Z[ak xk ](z) = Z[xk ](z/a). Calculamos Z[xk+1 ](z) = ∞ X xn+1 n=0 ∞ X n=0 ∞ X n=0 k0 −1 X n=0 xn z k0 −n . 22 ..
m ∈ N. se veriﬁca Z[k m xk ](z) = [−z d m ] Z[xk ](z). entonces desarrollando en serie de Laurent 1 1 1 = z−1 z1− 1 z X 1 1X 1 = = z n=0 z n n=0 z n+1 ∞ ∞ 100 . si xk = k 2 . entonces Z[k2 ](z) = [−z z d 2 d ] Z[1](z) = [−z ]2 dz µ 1 ¶ ¶ dz z − µ d d z d −z z2 = −z + −z = −z dz dz z − 1 dz z − 1 (z − 1)2 3z 2 2z 3 z − + . xn = Z −1 [F (z)]. donde 1 r ≥ 0. A. veamos que también lo es para m + 1. Por ejemplo. dz dz dz ∞ X d −xn nxn = z =z z n+1 dz z n n=1 n=1 ! ! Ã ∞ Ã∞ X xn d X xn d d = −z = z − x0 = −z Z[xk ](z). si F (z) = z−1 . Para calcular la transformada Z de una función F (z) basta calcular el desarrollo en serie de Laurent centrada en cero de manera que tenga un anillo de convergencia de la forma {z ∈ C : |z| > r}. entonces Z[kxk ](z) = ∞ X nxn n=0 ∞ X zn = Si suponemos el resultado cierto para m.2 Transformada Z inversa o equivalentemente Es interesante obtener transformadas Z inversas de funciones de variable compleja F (z). Para esto calculamos Z[k m+1 xk ](z) = Z[k · km xk ](z) = −z d Z[km xk ](z) dz d m d d = (−z )[−z ] Z[xk ](z) = [−z ]m+1 Z[xk ](z). Si m = 1. es decir. dz d donde por −z dz se entiende la operación derivada y luego multiplicación por −z. Hacemos la demostración por inducción en m. qué sucesiones veriﬁcan que Z[xn ](z) = F (z). Demostración. − n n dz z dz n=0 z dz n=1 ∞ X nxn n=1 zn Por ejemplo.Ecuaciones en diferencias • Dadas las sucesiones xk y km .2. = z − 1 (z − 1)2 (z − 1)3 si |z| > 1.
Por otra parte Z[1](z) = Entonces z .). Entonces la sucesión xk = Z −1 [1/(z − 1)] = (0.3 Aplicación a la resolución de la ecuación en diferencias lineales con coeﬁcientes constantes ½ yk+2 + yk+1 − 2yk = 1. y1 = 1. usando las propiedades de ésta y tomando en consideración las condiciones iniciales obtenemos Z[yk+2 + yk+1 − 2yk ](z) = Z[1](z). (z 2 + z − 2)(z − 1) −1 4 z2 3 = + . − 2 (z + 2) 2 (z − 1) (z − 1) z−1 z+2 y calculamos la transformada inversa obteniendo los desarrollos en series de Laurent ¶n X ∞ µ ∞ 1 1 1 X −2 (−2)n 1 = = = z + 2 z 1 − −2 z n=0 z z n+1 z n=0 si |z| > 2. 1. z−1 z−1 (z 2 + z − 2)Z[yk ](z) = z + con lo que Z[yk ](z) = Pasamos a fracciones simples Z[yk ](z) = z2 . A. Consideramos el problema obtenido anteriormente. . y desarrollando Z[yk+2 + yk+1 − 2yk ](z) = Z[yk+2 ](z) + Z[yk+1 ](z) − 2Z[yk ](z) = z 2 Z[yk ](z) − z + zZ[yk ](z) − 2Z[yk ](z) = (z 2 + z − 2)Z[yk ](z) − z. Finalmente ! Ã∞ ∞ ∞ X 1 X d 1 Xn+1 = =− z n+1 dz z n+1 z n+2 n=0 n=0 n=0 101 .. Tomando la transformada Z en la ecuación.2. z−1 z2 z = . 1. 1 1 1 = z−1 z1− −1 d = 2 (z − 1) dz µ 1 z−1 ¶ d = dz X 1 1X 1 = = z n=0 z n n=0 z n+1 ∞ ∞ 1 z si |z| > 1..Ecuaciones en diferencias si |z| > 1. y0 = 0. 1.
en que las raices son complejas: ½ xn+2 − 2xn+1 + 2xn = 1. = z −1 2z −1−i 2z −1+i 1 1 1 = z−1 z1− X 1 1X 1 = = . Para ello en primer lugar (z − 1)(z 2 z z = − 2z + 2) (z − 1)(z − 1 − i)(z − 1 + i) 1 1+i 1 1−i 1 − − . Si aplicamos la transformada Z a la ecuación. z−1 mientras que Z[xn+2 − 2xn+1 + 2xn ](z) = Z[xn+2 ](z) − 2Z[xn+1 ](z) + 2Z[xn ](z) = z 2 Z[xn ](z) − z 2 x0 − zx1 − 2zZ[xn ](z) − 2zx0 + 2Z[xn ](z) = (z 2 − 2z + 2)Z[xn ](z). de donde Z[xn ](z) = (z − 1)(z 2 z . x0 = x1 = 0. z n=0 z n n=1 z n ∞ ∞ Calculamos de forma separada 1 z 102 . z n=0 z n+2 ∞ X 4(−2)n+1 1 Xn+1 X 3 − − + z n=0 z n+2 z n+2 n=0 z n+2 n=0 por lo que si k ≥ 2 yk = 4(−2)k+1 − 4 + k. Entonces si |z| > 2 se tiene que Z[yk ](z) = −1 3 4 − + 2 (z − 1) z−1 z+2 ∞ ∞ ∞ Xn+1 X 1 X (−2)n −3 +4 = − z n+2 z n+1 z n+1 n=0 n=0 n=0 ∞ ∞ ∞ = 1 X −n − 4 + 4(−2)n+1 + = . Veamos a continuación el siguiente ejemplo. − 2z + 2) Desarrollamos la función en serie de Laurent para calcular xn .Ecuaciones en diferencias si |z| > 1. Por un lado ∞ X 1 1 = Z[1](z) = n z 1− n=0 1 z = z . tenemos que Z[xn+2 − 2xn+1 + 2xn ](z) = Z[1](z).
1+i z−1−i z1− z z n=0 z z n=1 con lo que agrupando ¶n X ∞ µ ∞ 1 1 1 1 1X 1−i = = = (1 − i)n−1 n . Por la fórmula integral de Cauchy 1 2πi Z γ Z[xk ](z)z n−1 1 dz = 2πi 1 2πi 103 (∗) = γ k=0 ∞ Z X k=0 γ Z X ∞ xk z n−k−1 dz xk z n−k−1 dz = xn . 1−i z−1+i z1− z z n=0 z z n=1 ∞ ∞ ∞ X 1 z 1 1X 1X n 1 − (1 + i) n − (1 − i)n n = 2 − 2z + 2) n (z − 1)(z z 2 n=1 z 2 n=1 z n=1 ¶ ∞ Xµ 1 1 1 n n 1 − (1 + i) − (1 − i) = . 4 4 nπ nπ − i sin ). . Supongamos una circunferencia del plano complejo γ que contiene todas las singularidades de la función Z[xk ](z)z n−1 .2. 2 2 zn n=1 y teniendo en cuenta que (1 + i)n = 2n/2 (cos (1 − i)n obtenemos nπ nπ + i sin ).4 Fórmula de inversión compleja Z[xk ](z) = ∞ X k=0 Supongamos que xk z −k . = 2n/2 (cos 4 4 ∞ y por tanto X³ z nπ ´ 1 n/2 = 1 − 2 cos .Ecuaciones en diferencias ¶n X ∞ µ ∞ 1 1 1 1 1X 1+i = = = (1 + i)n−1 n . de donde Z[xk ](z)z n−1 = ∞ X k=0 xk z n−k−1 . y multipliquemos ambos miembros de la igualdad por z n−1 . (z − 1)(z 2 − 2z + 2) n=1 4 zn xn = 1 − 2n/2 cos nπ . 4 A. para todo n ≥ 1.
y −3 que es de orden uno.Ecuaciones en diferencias dado que Z ½ xk z Z n−k−1 dz = γ 0 si n − k 6= 2. zi ) = ¢ dm−1 ¡ 1 lim m−1 (z − zi )m Z[xk ](z)z n−1 . i) = lim(z − i) z n−1 z→i (z − i)(z + i) n−1 z 1 1 = lim = in−1 = − in . (m − 1)! z→zi dz z(z − 1) . 2 z→−3 (z − 2) 25 lim (z − 3) De esta forma xn = Sea ahora 4 n 1 4 2 + n2n − (−3n ). = 25 25 50 y Res(f (z)z n−1 . Entonces Res(f (z)z n−1 . 25 50 25 f (z) = z2 1 . z→i (z + i) 2i 2 104 . m X i=1 Por el Teorema de los residuos 1 xn = 2πi γ Z[xk ](z)z n−1 dz = Res(Z[xk ](z)z n−1 . 2) = 1! z→2 dz (z − 2)2 (z + 3) µ n ¶ d z (z − 1) = lim z→2 dz (z + 3) n−1 (nz (z − 1) + z n )(z + 3) − z n (z − 1) = lim z→2 (z + 3)2 4 (n2n−1 + 2n )5 − 2n 1 = 2n + n2n . zi ) donde si zi es un polo de orden m. (z − 2)2 (z + 3) Vamos a modo de ejemplo a obtener la transformada inversa de f (z) = Sus polos son 2. −3) = z n (z − 1) z→−3 (z − 2)2 (z + 3) z n (z − 1) 4 = lim = − (−3n ). entonces el residuo es Res(Z[xk ](z)z n−1 . de orden dos. 2πixn si n − k = 2. Entonces µ ¶ d z n (z − 1) 1 n−1 2 lim (z − 2) Res(f (z)z . +1 cuyos polos simples son ±i.
obtenemos 1 1 xn = − in + (−i)n 2 2 π π −π −π n 1 1 + i sin ) = − (cos + i sin )n + (cos 2 2 2 2 2 2 1 nπ nπ 1 −πn −πn = − (cos + i sin ) + (cos + i sin ). 2 2 2 2 2 2 y dado que sin y cos obtenemos xn = − cos nπ πn = − sin 2 2 −πn nπ = − cos 2 2 nπ . 0 ≤ i ≤ n.2. Se deﬁne enotnces la función de transferencia asociada a la ecuación como T (z) = Z[yk ](z) 1 = . suponiendo que yi = 0 i < k. (A.5) siendo ai ∈ R. z→−i (z − i) 2i 2 lim (z + i) por lo que 1 1 xn = − in + (−i)n .Ecuaciones en diferencias y Res(f (z)z n−1 .. por lo que Z[yk ](z) = 1 an zn + an−1 z n−1 + .. 2 2 Expresando los números complejos en forma trigonométrica y utilizando la fórmula de De Moivre. Entonces. tomando la transformada Z obtenemos que (an z n + an−1 z n−1 + .5 Funciones de transferencia. La función de transferencia asociada a la transformada Z se deﬁne de forma análoga a la función de transferencia asociada a la transformada de Laplace. + a1 z + a0 Z[uk ](z). n+a n−1 + . Consideremos en este contexto una ecuación en diferencias ﬁnitas de la forma an yk+n + an−1 yk+n−1 + .... −i) = z n−1 z→−i (z − i)(z + i) n−1 z 1 1 = lim = − (−i)n−1 = (−i)n . + a1 yk+1 + a0 yk = xk .... + a1 z + a0 )Z[yk ](z) = Z[uk ](z). 2 A. + a z + a Z[uk ](z) an z n−1 z 1 0 105 .
zn+1 = f(n. es una ecuación autónoma y xn+1 = xn + yn . . zn . zn .. Con frecuencia. Esta ecuación puede reducirse a un sistema 1 2 k de orden uno de la manera siguiente.. Por ejemplo.. zn = yn+k−1 .. 106 ..5) es estable si y sólo si todos los polos de la función de transferencia veriﬁcan que |z| < 1. zn ) = (zn . Si la ecuación o sistema de orden uno no depende explícitamente de n.. . zn+1 )). zn+1 . yn .. entonces de forma compacta el sistema se puede escribir como donde f : Λ ⊆ Rk+1 → R es una función continua.. presentan lo que se conoce como comprotamiento caótico o complicado... . si zn = (zn .Ecuaciones en diferencias Podemos estudiar entonces la estabilidad de la ecuación entendiendo ésta de forma análoga al caso continuo estudiada en el tema anterior. zn ) donde 1 k 2 3 k 1 k f(n. Theorem 3 El sistema dado por la ecuación (A.. Las ecuaciones y sistemas de orden uno han tenido un reciente desarrollo debido a que son modelos de ciencias experimentales como la ecología y la economía.. ⎪ ⎩ k 1 k zn+1 = f (n.. es un sistema autónomo de orden uno. si para toda solución asociada a una condición inicial dada se veriﬁca que lim yk = 0.... zn = yn+1 . ⎨ .. ... zn+1 ).. yn+k−1 ).. es decir... f (n. zn .. k→∞ El siguiente resultado caracteriza la estabilidad del sistema en base a los polos de la función de transferencia.. .. yn+1 = 4yn (1 − yn−1 ). A...3 Ecuaciones no lineales yn+k = f (n. se dice que dicha ecuación o sistema es autónomo. yn+1 = xn − yn . . Entonces ⎧ 1 2 ⎪ zn+1 = zn . .. Deﬁnimos las variables zn = yn .. ⎪ k−1 k ⎪ z ⎪ n+1 = zn . Consideremos ahora una ecuación en diferencias de orden k no lineal 1 k Esto es. ⎪ 2 ⎪ 3 ⎪ zn+1 = zn .. zn+1 .. zn ).
Aguirre y P. V. S. Servicio de publicaciones. Myers. Braselton. Pascual. Paraninfo. Jiménez López. J. P. CRC Press.L. [Cas] [Dor] [Ela] [Jim] E. Faires. C. Ed. Wiley. Dormand. Fundamentos de ecuaciones diferenciales ( 2a edición). [NaSa] R. K.Bibliografía [AbBr] M. Numerical methods for ordinary diﬀerential equations. Abell y J. Cambridge University Press. [Wol] Stephen Wolfram. Ed. 1992. Elaydi. Thomson. [VJAP] L. Addison Wesley Longman. Ecuaciones diferenciales (con aplicaciones y notas históricas). Wolfram Media. [BuFa] R. J. Addison—Wesley Iberoamericana. curso universitario. M.P. Apuntes de ecuaciones diferenciales. McGraw—Hill. 2a Edición. [MaMy] B. N. [Oga2] K. Saﬀ. B. M. An introduction to diﬀerence equations. Métodos numéricos para la Física y la Ingeniería. [Sim] G. Química. D. Mathematical modeling of the Brusselator. R. S. 1995. Discrete—time control systems. D. [McPo] M. Métodos numéricos. ”The Mathematica Book”. 107 . McDowell y J. Numerical methods for diﬀerential equations. AP Proﬀessional. Prentice Hall. Universidad de Murcia. 2000. Jiménez. Nagle y E. Springer—Verlag. 1990. McGraw—Hill. Powers. F. Preprint Universidad de Notre Dame (2008). [Lam] J. Castillo y otros. Mathematica. Lambert. Mahan y R. J. Simmons. 1993. Diﬀerential Equations with Mathematica. Burden y J. Ogata. Vázquez.
Documents Similar To mastquim
Apuntes Sobre Ecuaciones Diferenciales
Ecuaciones_Diferenciales_-_UTP-2015-II_-1-__21147__
CAP1 Ecuaciones Diferenciales de PRIMER ORDEN Edwars Penney - 4Edi
ECUACION ORDINARIA.pdf
Ecuaciones diferenciales Economicas
Taller Uno 1-2-3 Ed Johnf Ecuaciones Diferenciales

References: resolución 
 resolución 
 resolución 
 resolución 
 Resolución 
 resolución