Esquema rk4 sistema

Esquema rk4 sistema gy DavidMuadDIb 1 110R6pp 16, 201 1 4 pagos M ‘todo de Runge-Kutta de cuarto orden e para un sistema de EDO Curso de F Isica Computacional M. en C. Gustavo Contreras May La aplicaci’ n del RK4 a un conjunto de EDO es an loga a la aplicaci ‘ n del m • todo de segundo o a o e orden.

Sea un conjunto de dos ecuaciones: y z El m» todo RK4 para ‘ste conjunto es e e kl = hf(yn , zn , tn ) II = hg(yn , zn , tn ) kl II h Q = hfyn+, zn * , +222k1 II h12 = hgyn hk3=hf 222 12 h k2 13 hg yn + zn + , tn +222k4—hf (yn *k3, zn + 13 , tn + h) 14 = hg(yn + k3, zn + 13 , tn + h) 1 yn+l yn + [kl +2k2 +2k3 + k4]61 zn+l = zn +[ll + 212 + 213 + 14]6 Ejercicio: Resuelve el siguiente conjunto de EDO con h = 0,3rty h – 0,51-1 y = z, z = -y, YO) = Izo -O -f ora to View nut*ge Explicaci’no Este programa se dise-o para resolver un conjunto de cualquier n mero de EDO de primer orden. u En el sub Swp to page subprograma FUNCT se definen el n’ mero de ecuaciones IM, asi como los IM valores de las u condiciones iniciales. Para utilizar el programa con un nuevo problema, el usuario debera cambiar las ecuaciones en FIJNCT, el valor de IM y las condiciones Iniciales. La estructura del programa es esencialmente la misma del programa RK4, pero se calcula cada paso intermedio en un ciclo DO para el n’ mero IM de ecuaciones. u 1 2. Variables Listado de variables Y(l) : y Y(2) : z Y(l) : I -‘sima inc gnita e o YN(I): yn para I = 1 y Zn para = 2, etc.

YA(I) : yn + kl 12 ‘ yn k2/ 2 yn + k3 /2 para 00 zn + II 12 zn + 12 12 zn + 13 12 para =200 II ,12 , 13, 14 1,2,3,4 : similar al anterior para la M sima ecuaci’ n diferencial e o IM : n» mero de ecuaciones en el conjunto u NS : n mero de intervalos de tiempo en un intervalo de impresi n, TD u XP : l’ mite m ‘ximo de t a TD : intervalo de impresi’ n para to C’ digo o A continuaci ‘ n Imite m’ ximo de t a TD : intervalo de impresi’n para to A continuaci’n se indica el c’ digo del programa. o DIMENSION PRINT * PRINT ‘Esquema RK4 para un conjunto de ecuaciones’ ! Numero de ecuaciones IM-2 ! Condicion inicial para yl en tzo Y(l) — 1 ! Condicion inicial 1 PRINT* PRINT *, ‘Intervalo de impresion para Y2 en tzo Y(2) = O de T? READ *, Pl PRINT *, ‘¿Numero de pasos en un intervalo de impresion de T? READ NS PRINT ‘T maximo para detener los aculos? READ XL H = PI,’NS PRINT ‘H – ‘ , HXP = O HH – H/2 PRINT * ! lnicializacion del numero de linea LI = O PRINT ‘Linea T YO) Y(2), ‘WRITE LI, xp, 28 LI = Ll+l DO N = 1 !

Tiempo anterior ! Tiempo nuevo ! Tiempo en el punto nuevo ! Esta parte calcula kl 2 DO 1, IM END DO XA=XB CALL ! Esta parte calcula Q DO I = 1, IM END DO XA=XM CALL ! Esta parte calcula k3 DO 3Lvf4 – 1, IM END DO XA-XM CALL ! ESta parte calcula kB DO 1=1, IM END DO XA=XM CALL ! Esta parte calcula k4 DO = 1, IM YA(I) = END DO XA=XP CALL DO I = 1, IM ! Esquerna RK4 Y(l) = Y(I) + END DO END DO WRITE («,98) LI, xp, MI), 98 FORMAT (IX, 12, FIO. 6, IP4E1 6. 6/ (15X, 1 P4El 6. )) IF (XP . LT. XL) GOTO 28 200 PRINT * PRINT *, ‘Oprime I para continuar o o para errnjnar READ KIF (K . EQ. l) GOTO I * END ! +++++++++++++++++++++++++++++++++++ + SUBROUTINE FUNCT(EK, J, YA. H) DIMENSION 0:10), ! Define un conjunto de ecuaciones EKO,I) = YA(2) * H EK(J, 2) = -YA(I) * H RETURN END N ‘tese que no est• implementado en el c’ digo, el almacenamiento de los datos en un archivo, por lo o a o que habr’ que agregarlo y posteriormente trabajar con ese archivo. a