EAIQ. Grado en Qu´ımica Pr´ actica 5. Curso 2014-15 Resoluci´ on num´ erica de EDOs con Matlab IMPORTANTE: Los alumnos deber´ an leer el texto de la pr´ actica antes de acudir a la clase correspondiente. 1. Introducci´ on Matlab dispone de un comando que permite resolver ecuaciones diferenciales, el comando dsolve ´ ticas II). En muchas ocasiones, este comando no resulta operativo ya que, o (ya visto en Matema bien no se puede encontrar la soluci´ on, o bien se obtiene una soluci´on impl´ıcita de la que no podemos extraer la informaci´ on que se buscaba. Un ejemplo de lo anterior es el problema de valor inicial: 0 √ y (x) = 0.1 y + 0.4 x2 , y(2) = 4. Si se intenta calcular una soluci´ on exacta haciendo uso de dicho comando: >> dsolve(’Dy=0.1*sqrt(y)+0.4*x^2’,’y(2)=4’,’x’) Warning: Possibly spurious solutions [solvelib::isolate] Warning: Explicit solution could not be found. > In dsolve at 101 ans = [ empty sym ] (La situaci´on no mejora si no ponemos la condici´on inicial). Supongamos que no sabemos resolver anal´ıticamente un problema de valor inicial con ninguno de los m´etodos estudiados, o bien que s´ı sabemos pero que resulta muy costoso. En cualquiera de estos casos, deberemos conformarnos con calcular una soluci´ on aproximada del problema y para ello resultan de gran utilidad los m´ etodos num´ ericos. Consideremos el problema de valor inicial dy = f (x, y), dx y(a) = ya (1) que supondremos tiene soluci´ on u ´nica, y(x), con x ∈ [a, b]. Un m´etodo num´erico nos va a permitir obtener un conjunto de puntos {xj , yj }nj=1 tales que yj ' y(xj ), xj ∈ [a, b], es decir, nos permite obtener el valor aproximado yj de la soluci´on buscada y(x) en una serie de puntos xj , j = 1, . . . , n, del intervalo [a, b]. El objetivo de esta pr´ actica es mostrar dos m´etodos num´ericos que permiten calcular una soluci´ on aproximada de un problema de valor inicial: el m´ etodo de Euler expl´ıcito y m´ etodos de tipo Runge-Kutta. Dichos m´etodos han sido estudiados en la clase te´orica, pero los recordamos aqu´ı. 1 2. M´ etodo de Euler Como acabamos de mencionar, nuestro objetivo es obtener una aproximaci´ on de la soluci´on del problema de valor inicial (1) en una serie de puntos (denominados nodos) del intervalo [a, b]. Uno de los m´etodos m´as sencillos para conseguir esto es el m´ etodo de Euler expl´ıcito. En primer lugar, determinaremos los puntos {xj } del intervalo [a, b] en los cuales vamos a aproximar la soluci´on. Para ello consideramos una partici´ on o discretizaci´ on del intervalo [a, b], es decir, una serie de puntos del intervalo que, por simplicidad, supondremos igualmente espaciados, verificando a = x1 < x2 < . . . < xn = b, de modo que xj = xj−1 + h, i = 2, . . . , n y h es un n´ umero menor que 1 que se denomina paso de discretizaci´ on y se obtiene como h= b−a n−1 (n es el n´ umero de nodos de discretizaci´ on). Cuanto m´as peque˜ no sea el paso, mayor ser´a el n´ umero de puntos en que aproximaremos la soluci´ on y, generalmente, mejor ser´a la aproximaci´on de la soluci´on. Si denotamos x1 , x2 , · · · , xn a los puntos de la partici´on, e y1 , y2 , · · · , yn a la aproximaci´on de la soluci´ on en dichos puntos, entonces y como se ha visto en clase, el m´ etodo de Euler calcula las aproximaciones de la siguiente manera: yj = yj−1 + f (xj−1 , yj−1 ) h, j = 2, · · · , n. (2) As´ı, y2 = y1 + f (x1 , y1 ) h siendo (x1 , y1 ) el punto sobre el que se da la condici´on inicial e y2 la aproximaci´on de la soluci´on en el punto x2 . A partir de ella calculamos y3 = y2 + f (x2 , y2 ) h obteniendo un nuevo par (x3 , y3 ) que representa el valor aproximado y3 de la soluci´on en el punto x3 , y as´ı sucesivamente. 2.1. Error de truncamiento o discretizaci´ on En la sucesi´on de valores y1 , y2 , . . . , yn generados a partir de la f´ormula (2), en general el valor de yj no coincidir´a con el de y(xj ) (salvo para j = 1), debido a que el algoritmo s´olo proporciona una aproximaci´on de la soluci´ on. As´ı pues, cometemos un error que mide la diferencia que hay entre la soluci´on exacta y la aproximada teniendo en cuenta los errores cometidos en cada paso del m´etodo. Se dice que el error es de orden α, lo que se denota como O(hα ), si existe una constante C independiente de h y un entero positivo α tales que |error|=|sol.exac - sol.aprox| ≤ Chα cuando h es suficientemente peque˜ na. As´ı, para el m´etodo de Euler el error es O(h). Esto implica que, si el tama˜ no del paso se reduce a la mitad, el nuevo error es, aproximadamente, C(h/2), es decir, el error se reduce, aproximadamente, en un factor 1/2. Esto podr´ıa llevarnos a pensar que si reducimos mucho el tama˜ no del paso, la aproximaci´ on ser´ a buen´ısima. Hay que tener cuidado pues esto no siempre es cierto ya que aqu´ı influyen mucho los denominados errores de redondeo, debidos, grosso modo, a que en el ordenador los n´ umeros s´ olo pueden representarse con un n´ umero finito de cifras significativas. Una manera de minimizar este error es reducir el n´ umero de operaciones a realizar y este n´ umero es tanto mayor cuanto m´ as peque˜ no sea h. Por otro lado, y dejando aparte los errores de redondeo, si h es demasiado peque˜ no el n´ umero de nodos aumenta considerablemente, lo que puede suponer un esfuerzo computacional enorme (a veces compensa una aproximaci´on peor pero obtenida en menor tiempo). 2 2.2. Ejemplos y ejercicios Ejemplo 1: La ley de enfriamiento de Newton establece que la raz´on de cambio de la temperatura T (t) de un cuerpo es proporcional a la diferencia entre la temperatura del cuerpo T (t) y la del entorno T∞ , esto es, 0 T (t) = k(T∞ − T (t)), donde k es la constante de transferencia de calor. Supongamos que k = 1 (min)−1 y que T∞ = 70o . Si el cuerpo est´a inicialmente a 100o , utilizar el m´etodo de Euler con un tama˜ no de paso de discretizaci´ on de 0.2 minutos para aproximar la temperatura del cuerpo al cabo de 0.8 minutos. Soluci´ on (Excel ): En Excel aproximaremos la soluci´on de un problema de valor inicial por el m´etodo de Euler como sigue. Ve´ amoslo para el ejemplo de la ley de enfriamiento de Newton. Introducimos en la columna A la partici´on del intervalo de tiempos (t = x), con paso 0,2 min. En la columna B, a la derecha del primer dato t1 = a = 0 min ponemos el correspondiente valor de y = T , o sea, la condici´ on inicial y1 = ya = Ta = 100o . Debajo de este valor hallamos y2 = T2 usando la f´ormula del m´etodo de Euler, es decir, tecleamos: = B1 + 1 ∗ (70 − B1) ∗ 0, 2 donde B1 es la celda que contiene a y1 = T1 , 1 = k, 70 = T∞ y h = 0, 2. Los siguientes valores T3 , ... se obtienen arrastrando hacia abajo a partir de la celda de T2 a la par de la columna de tj . Despu´es se representan gr´ aficamente. Observa la aproximaci´ on de T en t = 0,8 min. A continuaci´on la compararemos con la soluci´ on exacta, que la podemos obtener con Matlab como sigue: >> T=dsolve(’DT=70-T’,’T(0)=100’,’t’) % Calculamos la soluci´ on exacta T = 70+30*exp(-t) Introducimos esta funci´ on en Excel en la columna a la derecha de la soluci´on aproximada. Despu´es se representan ambas funciones simult´ aneamente seleccionando todos los datos y eligiendo insertar gr´afico en el men´ u. Ejercicio 1: Realiza el ejercicio anterior con paso de discretizaci´on de 0.1 minutos. Escribe a continuaci´on la estimaci´ on por Euler de T (1) y el valor exacto (a partir de la soluci´on exacta). Dibuja a continuaci´on las dos gr´ aficas. 3 Ejercicio 2: Consideremos el problema de valor inicial ( y = y 2 log(x); x y(1) = 1. y0 + (3) Calcular en el intervalo [1, 3] la soluci´ on exacta en Matlab con dsolve y luego las aproximaciones, en Excel , con el m´etodo de Euler con pasos 0.2, primero, y, luego, 0.1 en el intervalo [1, 3]. Representarlas en Excel . 3. M´ etodos de Runge-Kutta Uno de los m´etodos num´ericos m´ as eficientes para aproximar la soluci´on de un problema de valor inicial son los denominados m´ etodos de Runge-Kutta. Existen m´etodos de Runge-Kutta de distintos ´ordenes, pero nos centraremos aqu´ı en los de orden dos y cuatro. 3.1. M´ etodo de Runge-Kutta de segundo orden Aunque existen otras variantes en la familia, una forma de calcular las aproximaciones con el m´ etodo de Runge-Kutta de segundo orden es la siguiente: yj = yj−1 + h 2 f (xj−1 , yj−1 ) + f (xj , yj? ) , donde yj? = yj−1 + h f (xj−1 , yj−1 ) j = 2, · · · , n. 3.2. El comando ode23 Matlab dispone del comando ode23 que tiene programado una mejora del m´etodo de RungeKutta de segundo orden. As´ı, para resolver el problema de valor incial 0 y (x) = f (x, y), y(a) = ya , en el intervalo [a, b] mediante el comando ode23 tecleamos: [x,y]=ode23( f ,[a b] ,ya ); donde: f es el nombre de la funci´ on autom´ atica que contiene a la f (x, y) de la ecuaci´on diferencial; a y b son los extremos del intervalo donde se desea conocer la soluci´on; 4 ya es el valor de la soluci´ on en a (es decir el valor de la condici´on inicial y(a) = ya); x es un vector cuyas componentes contienen los nudos de discretizaci´on (los xi ); y es un vector cuyas componentes contienen los valores aproximados de la soluci´on en cada uno de los nodos xi (es decir, los valores yi ). Este comando no requiere como dato el paso de integraci´on h pues determina de manera autom´atica el tama˜ no del paso necesario para mantener los errores por debajo de una error predeterminado. En concreto, es 10−3 , para el error relativo, y 10−6 , para el error absoluto. Si se desea conocer la soluci´ on para ciertos valores de x, puede especificarse el paso de discretizaci´ on usando el comando con la opci´ on: [x,y]=ode23(f ,[a: h :b] ,ya) Ejemplo 2: Resolver el ejercicio 2 en Matlab pero usando la aproximaci´on de tipo Runge-Kutta dada por ode23. Soluci´ on Recordemos que en el ejercicio 2, se trataba de resolver el problema ( y y 0 + = y 2 log(x); x y(1) = 1. en el intervalo [1, 3]. Si utilizamos el comando ode23 para resolver con paso de discretizaci´on h = 0.2, las sentencias ser´ıan: >> f=@(x,y) y.^2.*log(x)-y./x >> [x23,y23]=ode23(f,[1:0.2:3],1) Ejercicio 3: Resuelve num´ericamente el ejercicio 2 utilizando el comando ode23 y con h=0.01. Representa gr´aficamente la soluci´ on exacta y la aproximaci´on. Anota a la izquierda las sentencias usadas y a la derecha esboza las dos gr´ aficas juntas y escribe el error m´aximo de la aproximaci´on. >> >> >> >> >> >> >> >> >> >> >> 5 4. EAIQ: Ejercicios pr´ actica ”Resoluci´ on num´ erica de EDOs”. Apellidos y nombre: . . . . . . . . . . . . . . . . . . Fecha y firma: . . . . . . . . . . . . . . . . . . . . . Grupo pr´ acticas: . . . . . . . Ejercicio 1: Se considera el problema de valor inicial siguiente: 0 √ y (x) = 0.1 y + 0.4x2 y(2) = 4 1. Aplica el m´etodo de Euler en Excel para hallar una soluci´on aproximada en [2,3] con h=0.2. Escribe la f´ormula tecleada y da el valor estimado en el punto x=3. 2. Repite el apartado anterior con h=0.1. 3. Esbozar las gr´ aficas, juntas, de las dos soluciones aproximadas. Ejercicio 2: Se supone que la temperatura de un cuerpo satisface la ecuaci´on diferencial dT = 2 · 10−6 (164 − T 4 ) dt y en el instante inicial del proceso es de 100o C. Se trata de estimar la temperatura al cabo de 5, 10 y 30 minutos. Hacerlo con el m´etodo de Euler en Excel con h=0.2, y con la variante de Runge-Kutta de 2o orden en Matlab , poniendo tambi´en h=0.2. Escribir, adem´as, la f´ormula tecleada en Excel y las sentencias usadas en Matlab . ¿Puedes calcular la soluci´ on exacta con dsolve? 6 Ejercicio 3: Un proceso de tercer orden A + 2B → productos tiene la ecuaci´on cin´etica dx = 0.05 (0.2 − x)(0.3 − 2x)2 , dt donde a = 0.2, b = 0.3 son cantidades iniciales del compuesto A y B, respectivamente, y k = 0.05 es la constante cin´etica. Teniendo en cuenta la condici´on inicial x(0) = 0, se pide : 1. Resolver num´ericamente la EDO en el intevalo [0, 100], utilizando los m´etodos de Euler y variante de Runge Kutta de 2o orden, para n = 100. Escribir la f´ormula de Excel , las sentencias usadas en Matlab y esbozar las gr´ aficas de las aproximaciones obtenidas. 2. Idem que el apartado anterior, variando la constante cin´etica (k = 0.5 y k = 0.005). ¿Puedes calcular la soluci´ on exacta con dsolve? 7
© Copyright 2024