Los métodos de Taylor de orden \(k\) tienen el inconveniente de tener expresiones más complejas a medida que va aumentando \(k\) debido a que las expresiones de la derivada \(k\)-ésima de la solución del problema de valores iniciales \(y^{(k)}\) es cada vez es más compleja y va aumentado de forma exponencial debido a la aplicación de la regla de la cadena.

Esto hace que los métodos de Taylor sean computancionalmente costosos para valores de \(k\) relativamente grandes y usualmente no se usen en la práctica.

Sin embargo, no hemos de olvidar que sabemos que el método de Taylor de orden \(k\) es consistente de orden \(k\), conocimiento que no siempre es tan sencillo de obtener en otro tipo de métodos.

En esta sección vamos a introducir los métodos de Runge-Kutta.

Veamos primero un ejemplo de método de Runge-Kutta introductorio: el método de Euler modificado: \[ \hat{y}_{i+1}=\hat{y}_i+h f\left(t_i+\frac{h}{2},\hat{y}_i+\frac{h}{2}f(t_i,\hat{y}_i)\right),\ \hat{y}_0 =y_0. \]

Observamos que no usamos las derivadas sucesivas de la función \(f(t,y(t))\) tal como hacíamos en los métodos de Taylor.

En este caso usamos el valor \(t_i+\frac{h}{2}\) y \(\hat{y}_i+\frac{h}{2}f(t_i,\hat{y}_i)\) como valores “cercanos” a \(t_i\) y \(\hat{y}_i\).

Calculemos el orden de consistencia del método anterior.

El error local de truncamiento será donde recordemos que \(y_i=y(t_i)\) es el valor de la solución exacta en el punto del mallado \(t_i\): \[ \begin{align*} \tau_{i+1} & =\frac{y_{i+1}-y_i}{h}-f\left(t_i+\frac{h}{2},y_i+\frac{h}{2}f(t_i,y_i)\right)\\ & =\frac{y_i+h y'(t_i)+\frac{h^2}{2}y''(t_i)+O(h^3)-y_i}{h}\\ &\quad -\left(f(t_i,y_i)+\frac{h}{2}\frac{\partial f}{\partial t}(t_i,y_i)+\frac{h}{2}\frac{\partial f}{\partial y}(t_i,y_i)f(t_i,y_i)+O(h^2)\right)\\ & = y'(t_i)+\frac{h}{2}y''(t_i)-f(t_i,y_i)-\frac{h}{2}\left(\frac{\partial f}{\partial t}(t_i,y_i)+\frac{\partial f}{\partial y}(t_i,y_i)f(t_i,y_i)\right)\\ & \quad +O(h^2) \end{align*} \]

Ahora, usando que \[ \begin{align*} y'(t_i) & =f(t_i,y_i),\\ y''(t_i) & =\frac{\partial f}{\partial t}(t_i,y_i)+\frac{\partial f}{\partial y}(t_i,y_i)y'(t_i) =\frac{\partial f}{\partial t}(t_i,y_i)+\frac{\partial f}{\partial y}(t_i,y_i)f(t_i,y_i), \end{align*} \] tenemos que \(\tau_{i+1}(h)=O(h^2)\).

Por tanto, el método de Euler modificado tiene orden de consistencia \(2\).

• INPUT valores inicial y final \(a\), \(b\), entero \(n\), condición inicial \(y_0\) y función \(f(t,y)\)
• Set \(h=\frac{b-a}{n}\).
• Set \(t=a\).
• Set \(\mathbf{\hat{y}}=\mathbf{0}\). (Definimos inicialmente el vector de aproximaciones como \(0\))
• Set \(\hat{y}_1=y_0\). (Definimos la primera componente del vector de aproximaciones como la condición inicial)

• For i=1,...,n
  • Set \(\hat{y}_{i+1} =\hat{y}_{i}+h f\left(t+\frac{h}{2},\hat{y}_{i}+\frac{h}{2}f(t,\hat{y}_{i})\right)\) (Hallamos el valor \(\hat{y}(t_{i})=\hat{y}(a+ih)\))
  • Set \(t=t+h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos el vector de aproximaciones).

Uno de los métodos más usados es el método de Runge-Kutta 4 cuya expresión es la siguiente: \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i +\frac{h}{6}(k_1+2k_2+2k_3+k_4),\\ k_1 & = f(t_i,\hat{y}_i),\\ k_2 & = f\left(t_i+\frac{h}{2},\hat{y}_i+\frac{h k_1}{2}\right),\\ k_3 & = f\left(t_i+\frac{h}{2},\hat{y}_i+\frac{h k_2}{2}\right),\\ k_4 & = f(t_i+h,\hat{y}_i + h k_3),\ \hat{y}_0=y_0 \end{align*} \] Observemos que sigue la filosofía de los métodos de Runge-Kutta, es decir, está definido en función de valores “cercanos” al valor “actual” \((t_i,\hat{y}_i)\) sin que intervengan las derivadas \(k\)-ésimas de la función \(f(t,y(t))\).

Proposición.

El orden de consistencia del método de Runge-Kutta 4 es \(4\).

La proposición anterior es la razón que el método se denomine método de Runge-Kutta 4.

• INPUT valores inicial y final \(a\), \(b\), entero \(n\), condición inicial \(y_0\) y función \(f(t,y)\)
• Set \(h=\frac{b-a}{n}\).
• Set \(t=a\).
• Set \(\mathbf{\hat{y}}=\mathbf{0}\). (Definimos inicialmente el vector de aproximaciones como \(0\))
• Set \(\hat{y}_1=y_0\). (Definimos la primera componente del vector de aproximaciones como la condición inicial)

• For i=1,...,n
  • Set \(k_1=f(t,\hat{y}_{i})\) (definimos los \(k_i, i=1,2,3,4\))
  • Set \(k_2=f\left(t+\frac{h}{2},\hat{y}_{i}+\frac{hk_1}{2}\right)\)
  • Set \(k_3=f\left(t+\frac{h}{2},\hat{y}_{i}+\frac{hk_2}{2}\right)\)
  • Set \(k_4=f\left(t+h,\hat{y}_{i}+h k_3\right)\)
  • Set \(\hat{y}_{i+1} =\hat{y}_{i}+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\) (Hallamos el valor \(\hat{y}(t_i)=\hat{y}(a+ih)\))
  • Set \(t=t+h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos el vector de aproximaciones).

Vamos a generalizar el método de Runge Kutta 4 visto anteriormente.

Un método de Runge-Kutta general es un método de un paso donde la expresión de la ecuación en diferencias tiene la expresión siguiente: \(\hat{y}_{i+1}=\hat{y}_i+h\varphi(t_i,\hat{y}_i),\) donde \[ \begin{align*} \varphi(t_i,\hat{y}_i) & =\sum_{j=1}^J c_j k_j,\\ k_j & = f\left(t_i+h a_j,\hat{y}_i+h\sum_{l=1}^J b_{jl} k_l\right),\ j=1,\ldots,J. \end{align*} \] También suponemos que \(\displaystyle\sum_{l=1}^J b_{jl}=a_j\) de cara a simplificar los cálculos a la hora de calcular el orden de consistencia del método.

Analicemos la expresión anterior:

• Los métodos de Runge-Kutta dependen de valores “cercanos” a \(t_i\): \(t_i+ha_j\), \(j=1,\ldots,J\) donde las \(a_j\) se han de determinar.
• Un método de Runge-Kutta se suele representar mediante el denominado tablero de Butcher:

• donde
  • en la primera columna aparece el vector \(\mathbf{a}=(a_1,\ldots,a_J)\) de coeficientes de \(h\) en la definición de las \(k_j\), \(j=1,\ldots,J\),
  • en la última fila aparece el vector \(\mathbf{c}=(c_1,\ldots,c_J)\) de coeficientes de \(k_j\) en la definición de \(\varphi(t_i,\hat{y}_i)\) y
  • en la tabla aparece la matriz \(\mathbf{B}=\begin{bmatrix}b_{11} & b_{12} & \ldots & b_{1J}\\ b_{21} & b_{22} & \ldots & b_{2J}\\ b_{J1} & b_{J2} & \ldots & b_{JJ}\end{bmatrix}\) de coeficientes de \(k_l\) en la definición de las \(k_j\), \(j=1,\ldots,J\).
• Si la matriz \(\mathbf{B}\) anterior de dimensiones \(J\times J\) es estrictamente triangular inferior con ceros en la diagonal, el método será explícito y, en caso contrario, será implícito.

• Recordemos que la suma de todas las filas excepto la primera componente nos da la primera columna: \(\displaystyle\sum_{l=1}^J b_{jl}=a_j\)

Vamos a escribir el método de Euler modificado según la notación anterior: \(\hat{y}_{i+1}=\hat{y}_i+h\varphi(t_i,\hat{y}_i),\) con \[ \begin{align*} \varphi(t_i,\hat{y}_i) & =k_2,\\ k_1 & = f(t_i,\hat{y}_i),\\ k_2 & = f\left(t_i+\frac{h}{2},\hat{y}_i+\frac{h}{2}k_1\right). \end{align*} \]

Por tanto, el valor de \(J\) vale \(J=2\) y el tablero de Butcher vale en este caso:

Como la matriz \(\mathbf{B}\) es triangular inferior con ceros en la diagonal, el método es explícito.

En el caso del método de Runge-Kutta 4, el valor de \(J\) vale \(J=4\) y el tablero de Butcher vale en este caso:

Como la matriz \(\mathbf{B}\) es triangular inferior con ceros en la diagonal, el método es explícito.

• Si uno quiere usar un método implícito de Runge-Kutta, el valor de \(J\) es pequeño debido a la complejidad del algoritmo.
• La ventaja de los métodos de Runge-Kutta es que no requieren el cálculo de las derivadas de la función \(f(t,y(t))\) tal como pasaba en los métodos de Taylor.
• La desventaja de dichos métodos es el cálculo del orden de la consistencia que suele ser bastante tedioso tal como hemos visto en el caso del método de Runge-Kutta 4.
• Para calcular el orden de la consistencia de un método de Runge-Kutta se tiene que aplicar el desarrollo de Taylor de una función de dos variables.

Si nos fijamos en los métodos de Euler modificado y de Runge-Kutta 4, la suma de las componentes del vector \(\mathbf{c}\) vale \(1\). La proposición siguiente nos dice que siempre suman \(1\) si queremos que el método sea consistente:

Consistencia en los métodos de Runge-Kutta.

Un método de Runge-Kutta es consistente si, y sólo si, \(\displaystyle\sum_{j=1}^J c_j=1\).

En el caso de los métodos de Runge-Kutta, la convergencia y la consistencia son equivalentes:

Proposición: convergencia en los métodos de Runge-Kutta.

Si un método de Runge-Kutta es convergente, es consistente.

Si un método de Runge-Kutta es consistente, explícito y además la función \(f\) del problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\) verifica la condición de Lipschitz respecto la segunda variable con constante de Lipschitz \(L\) para valores de \(h\) suficientemente pequeños, el método es convergente con el mismo orden que el orden de consistencia.

Hemos visto que un método de Runge-Kutta es consistente si, y sólo si, \(\displaystyle\sum_{j=1}^J c_j=1\).

Veamos qué condiciones aparecen si queremos consistencia de un determinado orden. Concretamente, veamos el resultado siguiente:

Proposición. Condiciones de consistencia para un método de Runge-Kutta.

Consideremos un método de Runge-Kutta \(\hat{y}_{i+1}=\hat{y}_i+h\varphi(t_i,\hat{y}_i),\) donde \[ \begin{align*} \varphi(t_i,\hat{y}_i) & =\sum_{j=1}^J c_j k_j,\\ k_j & = f\left(t_i+h a_j,\hat{y}_i+h\sum_{l=1}^{J} b_{jl} k_l\right),\ j=1,\ldots,J, \end{align*} \] con \(\displaystyle a_j=\sum_{l=1}^J b_{jl}\) para resolver el problema de valores iniciales \[ y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0. \]

Proposición. Condiciones de consistencia para un método de Runge-Kutta. (continuación)

Entonces,

1. si \(\displaystyle\sum_{j=1}^J c_j=1\), el método es consistente de orden como mínimo \(1\),
2. si además \(\displaystyle\sum_{j=1}^J a_j c_{j}=\frac{1}{2}\), el método tiene orden de consistencia como mínimo \(2\) y
3. si además \(\displaystyle\sum_{j=1}^J a_j^2 c_j = \frac{1}{3}\), \(\displaystyle\sum_{j,l=1}^J b_{jl}a_l c_j=\frac{1}{6}\), el método tiene orden de consistencia como mínimo \(3\).

Observaciones.

• La demostración anterior nos da un algoritmo para verificar órdenes de consistencia mayores. Sin embargo, la complejidad va aumentando a medida que el órden a verificar es mayor.

• La condición de orden de consistencia \(2\), \(\displaystyle \sum_{j=1}^J a_j c_j=\frac{1}{2}\), puede escribirse matricialmente \(\mathbf{c}^\top\mathbf{a}=\frac{1}{2}\), donde \(\mathbf{a}\) sería la primera columna del tablero de Butcher y \(\mathbf{c}\), la última fila de dicho tablero.

Observaciones. (continuación)

• Condiciones de orden de consistencia \(3\):
  • \(\displaystyle \sum_{j=1}^J a_j^2 c_j=\frac{1}{3}\) puede escribirse matricialmente como \(\mathbf{c}^\top\mathbf{a^2}=\frac{1}{3}\), donde \(\mathbf{a^2}=(a_1^2,\ldots,a_J^2)\) y el vector \(\mathbf{c}\) está indicado en la observación anterior.
  • \(\displaystyle\sum_{j,l=1}^J b_{jl}a_l c_j=\frac{1}{6}\), puede escribirse matricialmente \(\mathbf{c}^\top\cdot\mathbf{B}\cdot\mathbf{a}=\frac{1}{6}\) donde \(\mathbf{B}\) la matriz del tablero de Butcher y los vectores \(\mathbf{a}\) y \(\mathbf{c}\) están indicados en la observación anterior.

Todos los métodos numéricos para resolver el problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\) vistos hasta ahora se basan en introducir una ecuación en diferencias \(\hat{y}_{i+1}=\hat{y}_i +h\varphi(t,\hat{y}_i)\) con \(\hat{y}_0=y_0\) y hallar las aproximaciones \(\hat{y}_i =\hat{y}(t_i)=\hat{y}(a+h i)\), donde \(h\) es el paso definido por \(h=\frac{b-a}{n}\), siendo \(n\) el número de aproximaciones a realizar.

Dicho paso es fijo para todas las aproximaciones \(\hat{y}_i\), \(i=0,\ldots,n\).

En esta nueva sección vamos a considerar que el paso \(h\) puede cambiar para cada tiempo \(t_i\) dependiendo del comportamiento que tenga la solución aproximada \(\hat{y}(t_i)\) en dicho valor del tiempo \(t_i\).

La idea intuitiva es que si, fijado un paso \(h\), el error local de truncamiento en el tiempo \(t_i\), \(\tau_{i+1}(h)\), es grande, significa que hemos de avanzar con “pies de plomo” con pasos \(h\) pequeños; sin embargo, si el error local de truncamiento en el tiempo \(t_i\), \(\tau_{i+1}(h)\), es pequeño, podemos avanzar con pasos \(h\) relativamente grandes ya que tenemos margen para considerar dichos pasos.

Los algoritmos numéricos que tienen en cuenta el error local de truncamiento para calcular el paso \(h\) en cada tiempo \(t_i\) se denominan adaptativos y son los que introduciremos en esta sección.

Por tanto, no sabremos a priori los tiempos \(t_i\) donde vamos a calcular las aproximaciones \(\hat{y}_i\) sino que éstos se irán obteniendo dependiendo de los errores locales de truncamiento \(\tau_{i+1}(h)\) en el tiempo actual \(t_i\) calculando los pasos \(h\) “óptimos” en cada tiempo \(t_i\).

Concretamente, el algoritmo para calcular la aproximación final \(\hat{y}(b)\) en el tiempo final \(t=b\) sería el siguiente:

• Paso 0: consideramos \(t_i=a\).
• Paso 1: consideramos un valor \(h>0\).

• Paso 2: calculamos el valor real \(q>0\) tal que el error local de truncamiento en el paso \(qh\) sea menor que la tolerancia permitida \(\epsilon\), es decir, \(|\tau_{i+1}(qh)|\leq\epsilon\).
  • si \(q<1\), cambiamos \(h\) por \(qh\) y volvemos a empezar en el paso 1.
  • si \(q\geq 1\) y \(t_i+qh\leq b\), calculamos \(\hat{y}(t_i+qh)=\hat{y}(t_i)+qh\varphi(t_i,\hat{y}(t_i))\) y aceptamos como buena la aproximación \(\hat{y}(t_i+qh)\). Cambiamos el paso \(h\) por \(qh\), cambiamos el valor \(t_i\) por \(t_i + qh\) y volvemos al paso 1.
  • si \(q\geq 1\) y \(t+qh\geq b\), hacemos \(h=b-t\), calculamos \(\hat{y}(b)=\hat{y}(t)+h\varphi(t,\hat{y}(t))\) y acabamos hallando la aproximación \(\hat{y}(b)\).

Como podemos observar en el algoritmo anterior, el punto clave es estimar el error local de truncamiento \(|\tau_t(qh)|\) y ver si éste es o no menor que la tolerancia permitida \(\epsilon\).

Para estimar dicho error local de truncamiento consideramos dos métodos numéricos, \[ \begin{align*} \hat{y}_{i+1} & =\hat{y}_i +h\varphi(t_i,\hat{y}_i),\\ \tilde{y}_{i+1} & = \tilde{y}_i +h\tilde{\varphi}(t_i,\tilde{y}_i), \end{align*} \] el primero \(\hat{y}_i\), con un orden de consistencia \(k\) y el segundo, \(\tilde{y}_i\), con un orden de consistencia \(k+1\).

Por tanto, si \(y(t)\) es la solución exacta del problema de valores iniciales, podemos escribir: \[ \begin{align*} y(t_i+h) & =y(t_i)+h \varphi(t_i,y(t_i))+ O(h^{k+1}),\\ y(t_i+h) & =y(t_i)+h \tilde{\varphi}(t_i,y(t_i))+ O(h^{k+2}). \end{align*} \]

Entonces, si \(\tau_{i+1}(h)\) es el error local de truncamiento del primer método y considerando que \(y(t_i)\approx \hat{y}_i\), \[ \begin{align*} \tau_{i+1}(h) & =\frac{y(t_i+h)-y(t_i)}{h}-\varphi(t_i,y(t_i))=\frac{y(t_i+h)-\hat{y}_i}{h}-\varphi(t_i,\hat{y}_i)\\ & = \frac{y(t_{i+1})-(\hat{y}_i +h\varphi(t_i,\hat{y}_i))}{h}=\frac{1}{h}(y(t_{i+1})-\hat{y}_{i+1}). \end{align*} \] De la misma manera, si \(\tilde{\tau}_{i+1}(h)\) es el error local de truncamiento del segundo método y considerando que \(y(t_i)\approx \tilde{y}_i\), \[ \tilde{\tau}_{i+1}(h)=\frac{1}{h}(y(t_{i+1})-\tilde{y}_{i+1}) \]

Como consecuencia, \[ \begin{align*} \tau_{i+1}(h) & = \frac{1}{h}(y(t_{i+1})-\hat{y}_{i+1})=\frac{1}{h}((y(t_{i+1})-\tilde{y}_{i+1})+(\tilde{y}_{i+1}-\hat{y}_{i+1}))\\ & = \tilde{\tau}_{i+1}(h)+\frac{1}{h}(\tilde{y}_{i+1}-\hat{y}_{i+1}). \end{align*} \] Ahora bien, como \(\tau_{i+1}(h)=O(h^{k})\) y \(\tilde{\tau}_{i+1}(h)=O(h^{k+1})\) la “parte dominante” de la expresión anterior de \(\tau_{i+1}(h)\) viene de \(\frac{1}{h}(\tilde{y}_{i+1}-\hat{y}_{i+1})\) ya que el otro sumando tiene orden superior \(O(h^{k+1})\): \[ \tau_{i+1}(h)\approx \frac{1}{h}(\tilde{y}_{i+1}-\hat{y}_{i+1}). \] Fijémonos que la expresión anterior nos permite estimar el valor del error local de truncamiento del primer método ya que la parte de la derecha es computable a partir de los dos métodos considerados.

Como \(\tau_{i+1}(h)=O(h^k)\), existe una constante \(C\) tal que \(\tau_{i+1}(h)\approx C h^k\).

A continuación consideramos el paso \(qh\). Aplicando la expresión anterior que nos da el error local de truncamiento pero ahora para \(qh\), obtenemos: \[ \tau_{i+1}(qh)\approx C (qh)^k = q^k C h^k \approx q^k\tau_{i+1}(h)\approx \frac{q^k}{h}(\tilde{y}_{i+1}-\hat{y}_{i+1}). \] Recordemos que queremos calcular \(q\) tal que \(|\tau_{i+1}(qh)|\leq \epsilon\), siendo \(\epsilon\) la tolerancia permitida: \[ |\tau_{i+1}(qh)|\leq \epsilon,\ \Rightarrow \frac{q^k}{h}|\tilde{y}_{i+1}-\hat{y}_{i+1}|\leq \epsilon,\ \Rightarrow q\leq\left(\frac{\epsilon h} {|\tilde{y}_{i+1}-\hat{y}_{i+1}|}\right)^{\frac{1}{k}}. \]

La expresión anterior nos dice que

• si \(|\tilde{y}_{i+1}-\hat{y}_{i+1}|\) es pequeño, significa que tenemos margen para elegir el paso \(h\), ya que dos métodos numéricos, uno de orden superior a otro, nos da una diferencia pequeña en las aproximaciones. El valor \(q\) será relativamente grande,
• en cambio, si \(|\tilde{y}_{i+1}-\hat{y}_{i+1}|\) es grande, significa que en el tiempo \(t_i\) tenemos una diferencia grande en las aproximaciones de dos métodos numéricos de distinto orden lo que implica que el paso \(h\) elegido deber ser pequeño al haber “problemas” en las aproximaciones en zonas cercanas al tiempo \(t_i\). En este caso \(q\) será pequeño y hemos de reducir el paso \(h\) a \(qh\).

Para aplicar el método descrito anteriormente, tenemos que elegir dos métodos numéricos para resolver el problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\).

Si los métodos elegidos son un método de Runge-Kutta de orden de consistencia \(4\), \(O(h^4)\) y un método de Runge-Kutta con orden de consistencia \(5\), \(O(h^5)\), el método se denomina método de Runge-Kutta-Fehlberg.

Concretamente, las ecuaciones en diferencias de los métodos numéricos considerados son las siguientes: \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i +h\left(\frac{25}{216}k_1+\frac{1408}{2565}k_3+\frac{2197}{4104}k_4-\frac{1}{5}k_5\right),\\ \tilde{y}_{i+1} & = \hat{y}_i +h\left(\frac{16}{135}k_1+\frac{6656}{12825}k_3+\frac{28561}{56430}k_4-\frac{9}{50}k_5+\frac{2}{55}k_6\right), \end{align*} \]

con \[ \begin{align*} k_1 & = f(t_i,\hat{y}_i),\\ k_2 & = f\left(t_i+\frac{h}{4},\hat{y}_i +\frac{h}{4}k_1\right),\\ k_3 & = f\left(t_i+\frac{3h}{8},\hat{y}_i + h\left(\frac{3}{32}k_1+\frac{9}{32}k_2\right)\right),\\ k_4 & = f\left(t_i +\frac{12h}{13},\hat{y}_i+h\left(\frac{1932}{2197}k_1-\frac{7200}{2197}k_2+\frac{7296}{2197}k_3\right)\right),\\ k_5 & = f\left(t_i+h,\hat{y}_i +h\left( \frac{439}{216}k_1-8k_2+\frac{3680}{513}k_3-\frac{845}{4104}k_4\right)\right),\\ k_6 & = f\left(t_i+\frac{h}{2},\hat{y}_i+h\left( -\frac{8}{27}k_1+2k_2-\frac{3544}{2565}k_3+\frac{1859}{4104}k_4-\frac{11}{40}k_5\right)\right). \end{align*} \]

Observaciones.

• Las demostraciones que los dos métodos elegidos tienen órdenes de consistencia \(4\) y \(5\) se deja como ejercicio al lector. Son demostraciones tediosas que se tienen que realizar usando el procedimiento descrito en la sección donde se introducían los métodos de Runge-Kutta.
• Los métodos elegidos son métodos explícitos.
• Las expresiones \(k_i\) son las mismas para los dos métodos. Por tanto, dada la aproximación en el tiempo \(t_i\), \(\hat{y}_i\), las aproximaciones \(\hat{y}_{i+1}\) y \(\tilde{y}_{i+1}\) empiezan las dos con \(\hat{y}_i +\cdots\).
• Además, el tiempo de computación se reduce a la hora de calcular la estimación del error local de truncamiento que recordemos que vale: \[ \tau_{i+1}(h)\approx \frac{1}{h}(\tilde{y}_{i+1}-\hat{y}_{i+1}). \]

Observación.

Debido a que, para deducir el algoritmo, hemos hecho las aproximaciones \(y(t_i)\approx \hat{y}_i\) y \(y(t_i)\approx \tilde{y}_i\), se considera la siguiente expresión “ajustada y conservadora” del valor de \(q\): \[ q = \left(\frac{\epsilon h} {2|\tilde{y}_{i+1}-\hat{y}_{i+1}|}\right)^{\frac{1}{4}}\approx 0.84\left(\frac{\epsilon h} {|\tilde{y}_{i+1}-\hat{y}_{i+1}|}\right)^{\frac{1}{4}}= 0.84\left(\frac{\epsilon}{R}\right)^{\frac{1}{4}}. \]

• INPUT valores inicial y final \(a\), \(b\), paso mínimo \(h_\min\), paso máximo \(h_\max\), tolerancia permitida \(TOL\) condición inicial \(y_0\) y función \(f(t,y)\)
• Set \(t=a\). (Definimos el valor de \(t\) como el tiempo inicial \(a\))
• Set \(\hat{y}=y_0\). (Definimos el valor aproximado en el tiempo inicial como \(y_0\))
• Set \(h=h_\max\). (Consideramos el valor de \(h\) como el \(h\) máximo y lo iremos ajustando)

• While \(t\leq b\) (mientras el tiempo actual \(t\) no sobrepase el tiempo final, vamos avanzando)
  • Set \(k_1=f(t,\hat{y})\) (definimos los \(k_i, i=1,2,3,4,5,6\))
  • Set \(k_2=f\left(t+\frac{h}{4},\hat{y}+\frac{h}{4}k_1\right)\)
  • Set \(k_3=f\left(t+\frac{3h}{8},\hat{y}+ h\left(\frac{3}{32}k_1+\frac{9}{32}k_2\right)\right)\)
  • Set \(k_4=f\left(t +\frac{12h}{13},\hat{y}+h\left(\frac{1932}{2197}k_1-\frac{7200}{2197}k_2+\frac{7296}{2197}k_3\right)\right)\)
  • Set \(k_5=f\left(t+h,\hat{y}+h\left( \frac{439}{216}k_1-8k_2+\frac{3680}{513}k_3-\frac{845}{4104}k_4\right)\right)\)
  • Set \(k_6=f\left(t+\frac{h}{2},\hat{y}+h\left( -\frac{8}{27}k_1+2k_2-\frac{3544}{2565}k_3+\frac{1859}{4104}k_4-\frac{11}{40}k_5\right)\right)\)
  • Set \(R=\left|\frac{1}{360}k_1-\frac{128}{4275}k_3-\frac{2197}{75240}k_4+\frac{1}{50}k_5+\frac{2}{55}k_6\right|\) (Calculamos la estimación de \(|\tau_{i+1}(h)|\))

• • If \(R\leq TOL\) (si la estimación de \(|\tau_{i+1}(h)|\) es menor que la tolerancia, aceptamos la aproximación y avanzamos)
    • Set \(t=t+h\).
    • Set \(\hat{y}= \hat{y}+h\left(\frac{25}{216}k_1+\frac{1408}{2565}k_3+\frac{2197}{4104}k_4-\frac{1}{5}k_5\right)\).
  • Set \(q=0.84\cdot (TOL/R)^{\frac{1}{4}}\) (Vamos a ajustar el paso \(h\))

• • If \(q\leq 0.1\) (si \(q\) es muy pequeño, consideramos el paso \(h\) como \(0.1\cdot h\))
    • Set \(h=0.1\cdot h\).
  • Else
    • If \(q\geq 4\) (si \(q\) es muy grande, consideramos como \(h\), \(4h\) sin superar \(h_\max\))
      • Set \(h=\min\{4\cdot h,h_\max\}\).
    • Else
      • Set \(h=\min\{q\cdot h,h_\max\}\).

• • If \(t+h > b\) (Comprobamos si el “próximo” tiempo supera el tiempo final \(t=b\))
    • Set \(h=b-t\) (Si es el caso, ajustamos el paso \(h\) para llegar exactamente a \(t=b\))
  • Else (si el “próximo” tiempo no supera el tiempo final \(t=b\), miramos si es porque el paso \(h\) es demasiado pequeño)
    • If \(h<h_\min\)
      • Print Paso \(h\) demasiado pequeño. El algoritmo no funciona.
      • End (Acabamos).
    • If \(t=b\) (Comprobamos si hemos llegado al tiempo final \(t=b\))
      • Print \(\hat{y}\) (si es el caso damos el valor \(\hat{y}(b)\))