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