Subrutina RUKU

Ya se vió en que consiste el método de Runge-Kutta; ahora nos corresponde discutir como se realiza el cálculo en FORTRAN. La subrutina que realiza esta tarea es RUKU, con argumentos (Z, DZ), donde Z es el punto en que se inicia la integración y DZ corresponde a $y_{n+1} - y_{n}$; este no es precisamente el incremento sino la derivada. Podríamos entonces, escribir $DZ/h$ siendo entonces el incremento dado por hDZ, pero no tiene caso escribir h. DZ, es la derivada efectiva, no es realmente la derivada de la función en el punto Z; posiblemente sea más correcto hablar de una secante (en vez de una tangente a la curva) que multiplicada por el tiempo, nos da la función en un tiempo DT después, teniéndose con eso una aproximación a cuarto orden.

      SUBROUTINE  RUKU (Z,DZ)
      DIMENSION   Z(6),X(6),DZ(6),DX(6)
      COMMON      T,DT,H,ALFA,IPR,IPO,G1,G2,Z1,Z2,GPL,GMI,EPP,EPM,XK,MO
      COMMON      IC,IT,IG,IP
      CALL   ZERV (DZ)
      CALL   GRAD (Z,DX)
      CALL   AUGV (DZ,DZ,1.0/6.0,DX)
      CALL   AUGV (X,Z,DT/2.0,DX)
      CALL   GRAD (X,DX)
      CALL   AUGV (DZ,DZ,1.0/3.0,DX)
      CALL   AUGV (X,Z,DT/2.0,DX)
      CALL   GRAD (X,DX)
      CALL   AUGV (DZ,DZ,1.0/3.0,DX)
      CALL   AUGV (X,Z,DT,DX)
      CALL   GRAD (X,DX)
      CALL   AUGV (DZ,DZ,1.0/6.0,DX)
      RETURN
      END

Esta subrutina sí la vamos a discutir con todo detalle, tanto por su importancia como por estar formada con llamadas a otras subrutinas. Se comienza propiamente con la instrucción CALL ZERV(DZ) que realiza la operación DZ=O. Luego, con CALL GRAD(Z,DX) se está indicando la operación DX F(Z) $k_{0}$ que es la tangente a la curva en el punto $X_{0}$. No se escribe T en F porque para nuestro caso, la derivada no va a depender de esa variable. De la misma manera, vamos obteniendo CALL AUGV(DZ,DZ,1/6,DX) para obtener DZ=DZ+DX/6=DZ+$k_{0}/6$ CALL AUGV(X,Z,DT/2,DX) que es lo mismo que X=Z+(DT/2)*DX=Z+(DT/2)$k_0$ = $hF(y_n+k_0/2)$ Se ha estado tomando Y=Z, además no se está usando H en el cálculo de X. Ahora ya podemos calcular $k_{1}$: CALL GRAD(X,DX) o sea DX = F(X) = $k_{1}$ Así tenemos el siguiente paso: CALL AUGV(X,Z,1/3,DX) es decir X=Z+DX/3. Y ya tenemos dos términos del argumento de $k_{2}$. CALL AUGV(X,Z,DT/2,DX) equivalente a X = Z + DT/2(DX) = Z + (DT/2) $k_1$ Con eso ya podemos obtener $k_{2}$: CALL GRAD(X,DX) y tendremos así DX = F(X) $= k_{2}$ Para calcular $k_{3}$ que es la última $k$ que nos falta damos los siguientes pasos.

CALL AUGV(DZ,DZ,1/3,DX) para calcular DZ = DZ+DX/3 = DZ+ $k_{2}$/3 CALL AUGV(X,Z,DT,DX) X = Z+DT*DX = Z+DT$\cdot k_{2}$ CALL GRAD(X,DX) con esto tendremos $k_{3}$ DX = F(X) = F $( y_{n} + k_{2}) = k_{3}$ Y por último: CALL AUGV(DZ,DZ,1/6,DX) para obtener DZ=DZ+DX/6

Con lo cual queda completo un paso de la integración de Runge-Kutta en FORTRAN. Hay que observar que este método de Runge-Kutta a cuarto órden se reduce al de Simpson. Cuando hay dependencias del tiempo debe incluirse el argumento, calculando en términos de más múltiplos de DT.