´ ´ ´ PRACTICAS DE CALCULO NUMERICO AVANZADO ´ PRACTICA 3: M´ etodos etod os de Runge-Kutta Runge-Ku tta expl´ expl´ıcitos. ıcitos . Pares encajados enca jados En esta pr´actica actica implementaremos implementar emos un m´etodo etodo de Runge-Kutta Runge-Kut ta expl ex pl´´ıcito general que posteriormente post eriormente utilizare utilizaremos mos para construir construir una rutina adaptativ adaptativa a para integra integrarr sistemas sistemas de ecuacion ecuaciones es diferencia diferenciales les ordinarias (en concreto, implementaremos implementaremos el m´ etodo etodo de Runge-Kutta Fehlberg).
1
M´ etodos etod os de d e Runge-Kutta Run ge-Kutta expl´ expl´ıcitos
Construiremos un m´ etodo etodo de Runge-Kutta expl´ expl´ıcito general de s etapas etapas,, es decir, decir, que implem implemen ente te cualquier m´ etodo etodo de un paso con la estructura: s
yn+1 = yn + h
bi ki
(1)
i=1
donde
i−1
ki = f tn + ci h, yn + h
aij k j
, i = 1,...,s 1 ,...,s
(2)
j=1 j =1
Nuestros programas deber´an an ser v´alidos alidos para sistemas de ecuaciones diferenciales ordinarias de primer orden. Se sugiere sugiere utilizar utilizar alg´un un sistema de ecuaciones sencillo como comprobaci´on de funcionamiento de las rutinas que construyamos, en concreto, una posibilidad es:
y1 = y2 ; y2 = −y1
y1 (0) y2 (0)
=
0 1
(3)
Los programas deber´an an funcionar f uncionar para cualquier m´etodo etodo de Runge-Kutta sin m´as que proporcionar el tablero de Butcher correspondiente. Para facilitar la posterior tarea de construcci´on de un par encajado, ser´a conveniente realizar de forma separada el c´alculo de las etapas k i y el c´alculo alculo del paso (es decir, el c´ alculo alculo de yn+1 a partir de yn ). Construiremo Construiremoss la rutinas rutinas generales generales para m´ etodos etodos RK considerand considerando o como ejemplo ejemplo concreto concreto,, el siguiente m´ etodo etodo de orden 3:
c A bT
0 0 1/2 1/2 ≡ 3/4 0 2/9
0 0 3/4 1/3
0 0 0 4/9
(4)
En resumen, se pide construir las siguientes rutinas: 1. Funci´on on fsis: contenie conteniendo ndo la funci´on on f ( f (t, y ), donde y es el vector vector columna de las solucion soluciones. es. Por Por ejemplo, en la caso de la Ec. (3), podr p odr´´ıamos escribir: function f=fsis(t,y) f=[y(2);-y(1)];
2. Rutina Rutina etapas, en las que calcularemos las k i . En la misma rutin rutina a se deber´an an proporcionar las matrices c y A. La sintaxis ser´a: a: function k=etapas(tn,yn,h) k=etapas(tn,yn,h) %----------------------------------------------------------------
%Sintaxis: k=etapas(tn,yn,h) %---------------------------------------------------------------%funcion que calcula las etapas de un metodo de RK explicito para %sistemas de n EDOs de primer orden %Inputs: % tn: instante temporal % yn: vector solucion en el instante tn (vector columna n x 1) % h: paso %Outputs: % k: etapas (matriz n x s, donde s es el numero de etapas) %---------------------------------------------------------------%Dentro de la rutina se han de proporcionar las matrices c y A %---------------------------------------------------------------c=[0; 1/2; 3/4]; % Estas son las matrices c a=[0 0 0; % y A para un RK explicito de 3 1/2 0 0 ; % etapas y orden 3 0 3/4 0];
Las salidas de esta rutina son pues las ki definidas en la ecuaci´on (2). 3. Rutina paso: dadas los valores para las etapas k i , esta rutina implementa un paso del m´ etodo de Runge-Kutta correspondiente. En esta rutina se proporcionar´a el vector b del tablero de Butcher. La sintaxis y significado de inputs y outputs se describe a continuaci´on. function yout=paso(tn,yn,h) %---------------------------------------------------------------%Sintaxis: yout=paso(tn,yn,h) %---------------------------------------------------------------%funcion que calcula un paso para un metodo de RK del que se %proporciona el vector b y los valores de las etapas k_i %Inputs: % tn: instante temporal % yn: vector solucion en el instante tn (vector columna n x 1) % h: paso %Outputs: % yout: vector solucion en el instante temporal tn+h %---------------------------------------------------------------%Dentro de la rutina se ha de proporcionar la matriz b %---------------------------------------------------------------%Esta rutina llama a la funcion etapas(tn,yn,h) %---------------------------------------------------------------b=[2/9 1/3 4/9];%vector b para el RK de orden 3 dado como ejemplo
4. Rutina RK: esta rutina, utilizando rutina paso, integra la ecuaci´on diferencial definida por la funci´ on fsis con las condiciones iniciales y0. La sintaxis es la siguiente: function [t,y]=RK(t0,tf,y0,N) %--------------------------------------------------------%Sintaxis: [t,y]=RK(t0,tf,y0,N)
%--------------------------------------------------------%Integracion de un sistema de EDOs de primer orden %mediante un metodo de un paso %Inputs: % t0:instante inicial % tf:instante final % y0:vector columna con las n condiciones iniciales % N:numero de divisiones del intervalo [t0,tf] %Outputs: % t:vector conteniendo los instantes t0,t1,...,tN % y:matriz Nxn de las soluciones en t0,t1,...,tN %--------------------------------------------------------%Esta funcion debe llamar a una rutina para evaluar % y(n+1) a partir de y(n) mediante el correspondiente % metodo de un paso (por ejemplo, mediante un Runge-Kutta)
5. Rutina orden. Escribiremos un programa para comprobar que el orden del m´etodo de Runge Kutta que se da como prueba es, efectivamente, 3.
1.1
Algunos ejemplos de output como comprobaci´ on
Como ayuda para comprobar el funcionamiento de las rutinas, damos algunos ejemplos num´ericos para su comprobaci´on (aunque la mejor comprobaci´on ser´a la de obtener el orden mediante la rutina orden. >> format short e >> etapas(1,[0.1 0.9]’,0.01) ans = 9.0000e-01 -1.0000e-01
8.9950e-01 -1.0450e-01
8.9922e-01 -1.0675e-01
>> paso(1.4,[1 3]’,0.1) ans = 1.2945e+00 2.8852e+00 >> [t,y]=RK(0,1,[1 1]’,5) t = 0
2.0000e-01
4.0000e-01
6.0000e-01
8.0000e-01
1.0000e+00
1.0000e+00 1.0000e+00
1.1787e+00 7.8133e-01
1.3103e+00 5.3154e-01
1.3897e+00 2.6060e-01
1.4137e+00 -2.0704e-02
1.3813e+00 -3.0114e-01
y =
1.2
Listado de programas
En resumen, para esta primera parte se piden los siguientes programas: 1. orden.m 2. RK.m 3. etapas.m 4. paso.m 5. fsis.m
2
Programaci´ on de un par encajado
La idea de los pares encajados es simple: se basa en derivar pares de m´etodos de Runge-Kutta que compartan las mismas etapas ki y que sean de ´ordenes consecutivos. Es decir, supone disponer de un m´etodo de orden p s
yn+1 = yn + h
bi ki
i=1
y de un m´ etodo de orden p + 1:
s
yn+1 = yn + h
˜bi ki
i=1
donde
i−1
ki = f tn + ci h, yn + h
aij k j
, i = 1,...,s
j=1
ya que consideramos m´etodos expl´ıcitos. De este modo, se puede, por ejemplo, utilizar el m´ etodo de orden p para obtener la soluci´o n y el m´etodo de orden p + 1 para estimar el error local de truncamiento. De hecho, una estimaci´on del mismo para el m´ etodo de orden p est´ a dada por: ˜n+1 = h T
s
(bi − ˜bi )ki
(5)
i=1
Control adaptativo del paso de integraci´ on Si disponemos de una estimaci´on del error local de truncamiento, podemos plantearnos dise˜n ar un algoritmo que permita adaptar el paso de integraci´on. Si queremos asegurarnos de que ˜n+1 < Tol , T ˜n+1 siendo T ol la tolerancia al error local de truncamiento que queramos imponer, podemos estimar T como antes y comprobar si pasa o no el test. Si no lo pasa, reducimos el paso de integraci´on y volvemos a dar el paso. Pero, ¿c´omo ha de ser el nuevo paso que tomemos para proceder de forma eficiente?. Mediante argumentos simples podemos estimar que ese nuevo paso debe ser de la forma: hnuevo = q
T ol ˜n+1| |T
1/( p+1)
h
Hacer h=h ini
Calcular yn+1 ~ Calcular la estimacion Tn+1
(
Tol h n= q ~ Tn+1
n = n+1 ~ ¿Es Tn+1< Tol ?
1 p+1
)
hn
No
Si h n+1 =q
1 p+1
( TTol ) ~
hn
No
¿Es t n+1> b ?
n+1
Si Fin
Figura 1: diagrama de flujo para un par encajado siendo q, un factor de “seguridad” ( q ≈ 0.8 suele funcionar bien). En alg´ un momento lograremos un paso de integraci´on que supere el test. Cuando suceda esto, ser´ıa conveniente asegurarnos de que el siguiente paso que demos tambi´en va a resultar adecuado y econ´omico. Por tanto, avanzaremos el nuevo paso de integraci´on de acuerdo con la expresi´on anterior.
Algoritmo de un Runge-Kutta adaptativo Aunque los diagramas de flujo est´an algo pasados de moda, tienen su utilidad did´actica. Damos a continuaci´on el que corresponder´ıa a un m´etodo Runge-Kutta adaptativo t´ıpico:
Desarrollo de la pr´ actica El m´ etodo de Runge-Kutta-Fehlberg (4,5) responde al siguiente tablero de Butcher 0
0
1 4 3 8 12 13
1 4 3 32 1932 2197 439 216 8 − 27 25 216 16 135 1 360
1 1 2
bi ˜bi bi − ˜bi
0 0 9 32 − 7200 2197
−8 2 0 0 0
0 0 0 7296 2197 3680 513 − 3544 2565 1408 2565 6656 12825 128 − 4275
0 0 0 0 845 − 4104
1859 4104 2197 4104 28561 56430 2197 − 75240
0 0 0 0 0 − 11 40 − 15 9 − 50 1 50
0 0 0 0 0 0 0 2 55 2 55
Este m´ etodo utiliza un m´etodo de cuarto orden de 5 etapas y un m´ etodo de quinto orden de 6 etapas. Se proporcionan los coeficientes en formato Matlab en: http://personales.unican.es/segurajj/coefs.txt
Implementaci´ on del m´ etodo Comenzaremos por repetir la primera parte de la pr´actica para los m´ etodos de Runge-Kutta de orden 4 y 5 descritos anteriormente, comprobando el orden de convergencia. Para ello, modificaremos ligeramente las rutinas etapas, y paso.
1. etapas45: igual que etapas.m pero con los coeficientes A correspondientes al m´etodo de RungeKutta Fehlberg. La sintaxis ser´a la misma, salvo porque, para aplicar f´acilmente las rutinas a varias funciones, an˜adiremos una primera entrada para especificar la funci´on f (t, y) correspondiente. La sintaxis ser´a entonces: function k=etapas45(fun,t,yn,h)
Siendo fun una cadena alfanum´erica. Por ejemplo, fun=’fsis’, si queremos considerar la f (t, y) definida en el fichero fsis.m. Otra diferencia respecto a etapas.m, ser´a en la propia evaluaci´on de las etapas para la que se utilizar´ a el comando feval. Recordemos que, por ejemplo, feval(fun,1,[3 -1]’) para fun=’fsis’ equivale a escribir fsis(1,[3 -1]’) . Para m´as detalles, consultar el “help”. 2. paso45: igual que paso.m pero incluyendo los coeficientes b para los Runge-Kutta encajados. Las salidas de la rutina ser´an el nuevo valor dado por el m´ etodo de orden 5 as´ı como la estimaci´on del error (ec. (5)): function [yout,err]=paso45(fun,tn,yn,h) %------------------------------------------------------------%Sintaxis: [yout,err]=paso45(fun,tn,yn,h) %-------------------------------------------------------------%funcion que calcula un paso para el metodo de RK-F del que %se proporcionan los vectores b y los valores de las etapas k_i %Inputs: % fun: cadena alfanumerica con el nombre del fichero que % contiene la funcion f(t,y) % tn: instante temporal % yn: vector solucion en el instante tn (vector columna n x 1) % h: paso %Outputs: % yout: vector solucion en el instante temporal tn+h utilizando % el RK de orden 5 % err: estimacion del error de truncamiento local para el RK-F %---------------------------------------------------------------%Dentro de la rutina se han de proporcionar las matrices b (para %orden (4) y bg (para orden 5)) %---------------------------------------------------------------%Esta rutina llama a la funcion etapas45(tn,yn,h) %---------------------------------------------------------------b=[25/216 0 1408/2565 2197/4104 -1/5 0]; bg=[16/135 0 6656/12825 28561/56430 -9/50 2/55];
3. RK45, igual que RK.m pero para el Runge-Kutta Fehlberg. El vector de salida ser´a la soluci´on num´ erica para el m´ etodo de orden 5. Sintaxis: [t,y,error]=RK45(fun,t0,tf,y_0,N)
donde la u ´ nica diferencia respecto a RK.m es la aparici´on del primer input, antes descrito. 4. orden45.m: para comprobar que el orden obtenido con RK45 es, efectivamente, 5.
Una vez verificado el correcto funcionamiento de RK45, utilizaremos las funciones etapas45 y etodo adaptativo de Runge-Kutta Fehlberg, que implementa el diagrama de paso45 para construir el m´ flujo de la Figura 1. La sintaxis ser´a la siguiente: function [t,y,errlt]=fehlb(fun,t0,tf,y_0,eps) %-------------------------------------------------------------------%Sintaxis: [t,y,errlt]=fehlb(fun,t0,tf,y_0,eps) %-------------------------------------------------------------------%funcion que implementa el metodo de Runge-Kutta Fehlberg adaptativo %Inputs: % fun: cadena alfanumerica para el nombre del fichero conteniendo la % funcion f(t,y) % tn: instante inicial % tf: instante final % y_0: vector columna de condiciones iniciales % eps: tolerancia de error para el error local de truncamiento %Outputs: % t: vector conteniendo los instantes temporales % y: valores numericos proporcionados para los instantes t; es una % matriz n x N (n: numero de ecuaciones; N: longitud del vector t) % err: estimacion del error de truncamiento local para cada instante %--------------------------------------------------------------------%Dentro de la rutina se han de proporcionar las matrices b (para orden %4) y bg (para orden 5)) %--------------------------------------------------------------------%Esta rutina llama a la funcion paso45(fun,tn,yn,h) %---------------------------------------------------------------------
Para la construcci´o n del la rutina fehlb.m se puede tomar cualquier paso inicial (menor que la longitud del intervalo) puesto que la rutina adaptar´a el h hasta que resulte conveniente para la precisi´on exigida. Por ejemplo, tomaremos h = (tf − t0)/10. Adem´ as de las rutinas etapas45, paso45, RK45 y fehlb, construiremos tres funciones que contengan distintas f (t, y): 1. fsis: que hemos dados expl´ıcitamente en la primera p´agina de este gui´on. 2. fpen: que implementa la f (t, y) para el p´ endulo simple (tomemos g = 9.8ms
2,
−
L = 1m).
3. fbru: que implementa la f (t, y) para el sistema conocido como Brusselator: y1 = 1 + y12 y2 − 4y1 y2 = 3y1 − y12 y2
(6)
sistema para el que, a diferencia de los otros sistemas estudiados, no tenemos soluci´on anal´ıtica, ni siquiera en t´ erminos de funciones especiales como ocurre con el p´ endulo simple. Adem´ as, para poder comparar en el caso del p´ endulo, construiremos la rutina analitica.m , que proporciona la soluci´on anal´ıtica del p´ endulo simple considerada en la anterior pr´actica. La sintaxis ser´a: function theta=analitica(t,cond_in);
Esta rutina se proporciona en http://personales.unican.es/segurajj/analitica.txt A continuaci´on, describiremos las actividades a realizar con los programas desarrollados en la segunda parte de esta pr´actica. Esta actividades deber´an ser ejecutadas por un programa Matlab al que llamaremos adaptativo.m .
2.1
Listado de programas
En resumen, los programas necesarios para realizar esta parte de la pr´actica son: 1. orden45.m 2. RK45.m 3. etapas45.m 4. paso45.m 5. adaptativo.m 6. fehlb.m 7. fsis.m (programa ya utilizado en la primera parte) 8. fpen.m 9. fbrus.m 10. analitica.m (programa que se proporciona en la p´agina web).
2.2 2.2.1
Actividades a realizar (adaptativo.m) Primer problema
Resolveremos el problema de valores iniciales de le Ec. (3) para t ∈ [0, 10]. Dibujaremos cuatro gr´aficas (que se recomienda que aparezcan en la misma ventana gr´afica utilizando subplot): 1. Aplicaremos el m´ etodo adaptativo para una tolerancia de error absoluto = 10 10 . Dibujaremos en una misma gr´afica la soluci´o n y el paso h como funci´on del tiempo. Como los pasos ser´an peque˜ nos, se recomienda multiplicarlos por un factor para que sean visibles en la gr´afica (un factor 20 debe bastar). −
2. Dibujaremos el error absoluto cometido para la primera componente del vector soluci´on, comparando con la soluci´on exacta y1 (t) = sin t. 3. Considerando el mismo n´umero de puntos que resultan en el m´ etodo adaptativo pero escogiendo un paso constante, resolveremos el mismo problema mediante RK45 y dibujaremos la soluci´on. 4. Representar el error para la primera componente obtenida mediante RK45 comparando con la soluci´ on exacta. ¿Tiene en este caso alguna ventaja el m´etodo adaptativo RK-Fehlberg?, ¿Por qu´e?.
2.2.2
Segundo problema (p´ endulo: primera parte)
Resolveremos el problema del p´ endulo simple compar´andolo con la soluci´on anal´ıtica. Las condiciones iniciales ser´a n: ´angulo inicial π/2 y velocidad angular inicial nula. Consideraremos el intervalo t ∈ [0, 10] y una tolerancia para el error absoluto = 10 10 . Repetimos las representaciones gr´aficas consideradas para el primer problema. La soluci´on anal´ıtica en este caso se obtiene mediante analitica.m . Los pasos h convendr´a escalarlos por un factor ∼ 100 para que aperezcan claramente en la gr´afica. Adem´ as, en otra ventana gr´afica proporcionaremos dos gr´aficas con las estimaciones de los errores de truncamiento locales para los dos m´ etodos ( fehlb y RK45) −
2.2.3
Tercer problema (p´ endulo: segunda parte)
Repetimos las gr´aficas del segundo problema pero tomamos como condiciones iniciales: angulo inicial 0.999999π, velocidad angular inicial 0. Consideraremos el intervalo t ∈ [0, 20] y una tolerancia para el error absoluto = 10 14 . Al igual que antes, utilizaremos un factor ∼ 100 para los pasos h. ¿Qu´ e diferencias se observan respecto a los anteriores valores iniciales?. −
2.2.4
Cuarto problema: el Brusselator
Consideraremos el problema del Brusselator con condiciones iniciales y 1 (0) = 1.5, y2 (0) = 3.1 en el intervalo t ∈ [0, 30]. Como tolerancia de error absoluto tomaremos = 10 10 . Dibujaremos dos gr´aficas en una misma ventana gr´afica: −
1. Gr´afica con las soluciones y1 , y2 y los pasos h en funci´on de t (se recomienda multiplicar los pasos por ∼ 50 para que sean visibles en la gr´afica). 2. Gr´afica con la estimaci´on del error local de truncamiento. Discutir brevemente los resultados.
