Los métodos vistos hasta ahora son métodos de un paso debido a que para calcular la solución aproximada \(\hat{y}_{i+1}=\hat{y}(t_{i+1})\) en \(t=t_{i+1}\) sólo usamos el valor de la solución aproximada en el tiempo anterior \(t_i\), \(\hat{y}_i=\hat{y}(t_i)\).

Ahora bien, como en general, el error cometido entre los valores de la solución aproximada \(\hat{y}_i=\hat{y}(t_i)\) y los valores de la solución exacta \(y_i=y(t_i)\), \(|y_i-\hat{y}_i|\) tienden a crecer a medida que el índice \(i\) aumenta, o a medida que nos acercamos al tiempo final \(t=b\), parece razonable de cara a que dicho error no crezca demasiado desarrollar métodos numéricos que, de cara a calcular el valor de la solución aproximada \(\hat{y}_{i+1}=\hat{y}(t_{i+1})\) usen no sólo el valor de la solución aproximada en el tiempo anterior \(t_i\), \(\hat{y}_i=\hat{y}(t_i)\) sino valores previos a éste, \(\hat{y}_{i-1}=\hat{y}(t_{i-1}),\hat{y}_{i-2}=\hat{y}(t_{i-2})\ldots,\hat{y}_{i-m}=\hat{y}(t_{i-m})\) para un \(m>0\) fijado.

Dichos métodos se denominan métodos multipaso cuya definición formal es la siguiente:

Consideremos el problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\).

Sea \(n\) un entero positivo y \(h=\frac{b-a}{n}\) el paso elegido. Consideramos el mallado \(t_i=a+i\cdot h\), \(i=0,\ldots,n\). Diremos que un método numérico es un método multipaso de orden \(m\) si la ecuación en diferencias correspondiente es de la forma: \[ \begin{align*} \hat{y}_{i+1} & = a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m}+\\ &\quad h(b_m f(t_{i+1},\hat{y}_{i+1})+b_{m-1}f(t_i,\hat{y}_i)+\cdots +b_0 f(t_{i+1-m},\hat{y}_{i+1-m})), \end{align*} \] para \(i=m-1,m,\ldots,n-1\), donde se tienen que especificar los valores iniciales siguientes: \[ \hat{y}_0=y_0,\ \hat{y}_1=\tilde{y}_1,\ldots,\hat{y}_{m-1}=\tilde{y}_{m-1}. \]

Observaciones.

• Un método multipaso tiene dos partes,
  • la parte lineal definida a partir de los coeficientes \(a_i\), \(i=0,\ldots,m-1\) y
  • la parte no lineal definida a partir de los coeficientes \(b_i\), \(i=0,\ldots,m-1\).
• Si \(b_m=0\), el método multipaso es explícito y, en caso contrario implícito. En este segundo caso, el valor \(\hat{y}_{i+1}\) aparece en los dos miembros de la ecuación y se tiene que resolver una ecuación no lineal para hallar dicho valor para cada valor de \(i\). Esto hace que la complejidad aumente. Sin embargo, los métodos implícitos suelen tener un orden de consistencia mayor que los explícitos.

Observaciones. (continuación)

• Para poder aplicar un método multipaso, necesitamos \(m\) valores iniciales, \(\hat{y}_0,\hat{y}_1,\ldots,\hat{y}_{m-1}\). A partir de dichos valores, vamos hallando los siguientes valores aproximados: \[ \begin{align*} \hat{y}_{m} & = a_{m-1}\hat{y}_{m-1} +a_{m-2}\hat{y}_{m-2}+\cdots +a_0 \hat{y}_{0}+\\ &\quad h(b_m f(t_{m},\hat{y}_{m})+b_{m-1}f(t_{m-1},\hat{y}_{m-1})+\cdots +b_0 f(t_{0},\hat{y}_{0})),\\ \hat{y}_{m+1} & = a_{m-1}\hat{y}_{m} +a_{m-2}\hat{y}_{m-1}+\cdots +a_0 \hat{y}_{1}+\\ &\quad h(b_m f(t_{m+1},\hat{y}_{m+1})+b_{m-1}f(t_{m},\hat{y}_{m})+\cdots +b_0 f(t_{1},\hat{y}_{1})),\\ &\vdots\\ \hat{y}_{n} & = a_{m-1}\hat{y}_{n-1} +a_{m-2}\hat{y}_{n-2}+\cdots +a_0 \hat{y}_{n-m}+\\ &\quad h(b_m f(t_{n},\hat{y}_{n})+b_{m-1}f(t_{n-1},\hat{y}_{n-1})+\cdots +b_0 f(t_{n-m},\hat{y}_{n-m})). \end{align*} \]

Observaciones. (continuación)

• Para calcular los valores iniciales \[ \hat{y}_0=y_0,\ \hat{y}_1=\tilde{y}_1,\ldots,\hat{y}_{m-1}=\tilde{y}_{m-1}, \] se suele usar un método de un paso como por ejemplo Runge-Kutta 4.

• Para ver que un método multipaso tiene orden de consistencia \(k\), hay que demostrar la expresión siguiente tiene \(O(h^k)\): \[ \begin{align*} \tau_{i+1}(h) & =\frac{y_{i+1}-(a_{m-1}y_i +a_{m-2}y_{i-1}+\cdots +a_0 y_{i+1-m})}{h}\\ &\quad -(b_m f(t_{i+1},y_{i+1})+b_{m-1}f(t_i,y_i)+\cdots +b_0 f(t_{i+1-m},y_{i+1-m})\\ & =O(h^k), \end{align*} \] donde recordemos que \(y_i=y(t_i)\), el valor de la solución exacta del problema de valores iniciales en \(t=t_i\).

Observaciones. (continuación)

• Seguidamente tenemos que desarrollar por Taylor las funciones \(y_{i+1},y_{i-1},\ldots,y_{i+1-m}\) alrededor de \(t=t_i\) y las funciones de dos variables \(f(t_{i+1},y_{i+1}),f(t_{i-1},y_{i-1}),\ldots,f(t_{i+1-m},y_{i+1-m})\) alrededor del punto \((t_i,y_i)\) y ver que coeficientes de las potencias \(h^j\) se anulan. Entonces, si los coeficientes de \(h^0, h^1,\ldots,h^{k-1}\) son cero, el orden de consistencia del método será \(k\).

Vamos a ver una metodología para deducir métodos multipaso basada en usar fórmulas de interpolación numérica para integrar una función real de variable real en nodos equiespaciados. Ver el curso Métodos numéricos con Python. Cálculo numérico para más detalles.

Recordemos el problema de valores iniciales a resolver: hallar una función \(y(t)\) tal que \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\).

Si consideramos el intervalo \([t_i,t_{i+1}]\) donde \(t_i=a+i\cdot h\) y \(t_{i+1}=a+(i+1)\cdot h\), \(h=\frac{b-a}{n}\) son dos puntos del mallado, podemos escribir: \[ y_{i+1}-y_i =y(t_{i+1})-y(t_i)=\int_{t_i}^{t_{i+1}}y'(t)\, dt. \] Ahora bien, como \(y'(t)=f(t,y(t))\) la expresión anterior será: \[ y_{i+1}-y_i = \int_{t_i}^{t_{i+1}}f(t,y(t))\, dt. \]

A continuación, consideramos el polinomio interpolador de grado \(m\) de la función \(f(t,y(t))\), \(P_m(t)\), en los nodos equiespaciados del mallado \(t_i, t_{i-1},\ldots t_{i+1-m}\). Podemos escribir la expresión siguiente: \[ f(t,y(t))=P_m(t)+\frac{f^{(m)}(\xi_i(t),y(\xi_i(t)))}{m!}(t-t_i)(t-t_{i-1})\cdots (t-t_{i+1-m}), \] con \(\xi_i(t)\in (t_i,t)\).

Como los nodos \(t_j\), \(j=i,i-1,\ldots,i+1-m\) son equiespaciados, tenemos que el polinomio interpolador se puede escribir de la forma siguiente basado en las diferencias hacia atrás: \[ P_m(t)=\sum_{k=0}^{m-1}(-1)^{k}\binom{s}{k}\nabla^{k} f(t_i,y(t_i)), \]

donde \(s\) es la nueva variable que vale \(s=\frac{t_i-t}{h}\), \(\binom{s}{k}=\frac{s(s-1)\cdots (s-k+1)}{k!}\), y \[ \begin{align*} \nabla^0 f(t_i,y(t_i)) & =f(t_i,y(t_i)),\\ \nabla^k f(t_i,y(t_i)) & =\nabla^{k-1}f(t_{i},y(t_i))-\nabla^{k-1}f(t_{i-1},y(t_{i-1})). \end{align*} \]

Fijémonos que, como \(t\in [t_i,t_{i+1}]\), entonces \(s\in [-1,0]\). Por tanto, haciendo el cambio de variable \(s=\frac{t_i-t}{h}\), \(dt=-hds\), tenemos que: \[ \begin{align*} y_{i+1}-y_i & = \int_{t_i}^{t_{i+1}}f(t,y(t))\, dt =-h\int_0^{-1}\sum_{k=0}^{m-1}(-1)^{k}\binom{s}{k}\nabla^{k} f(t_i,y(t_i))\,ds\\ & \quad +\frac{h^{m+1}(-1)^{m+1}}{m!}\int_{0}^{-1}f^{(m)}(\xi_i(s),y(\xi_i(s)))s (s-1)\cdots (s-(m-1))\, ds, \end{align*} \] con \(\xi_i(s)\in (-1,0)\).

En las integrales anteriores, hacemos otro cambio de variable \(\hat{s}=-s\), con \(ds=-d\hat{s}\) y nos queda: \[ \begin{align*} y_{i+1}-y_i & = h\int_0^{1}\sum_{k=0}^{m-1}(-1)^{k}\binom{-\hat{s}}{k}\nabla^{k} f(t_i,y(t_i))\,d\hat{s}\\ & \quad +\frac{h^{m+1}}{m!}\int_{0}^{1}f^{(m)}(\hat{\xi}_i(\hat{s}),y(\hat{\xi}_i(\hat{s}))\hat{s} (\hat{s}+1)\cdots (\hat{s}+(m-1))\, d\hat{s} \\ & = h \sum_{k=0}^{m-1} \nabla^{k} f(t_i,y(t_i))(-1)^k \int_0^1 \binom{-\hat{s}}{k}\,d\hat{s}\\ & \quad +\frac{h^{m+1}}{m!}\int_{0}^{1}f^{(m)}(\hat{\xi}_i(\hat{s}),y(\hat{\xi}_i(\hat{s})))\hat{s} (\hat{s}+1)\cdots (\hat{s}+(m-1))\, d\hat{s}, \end{align*} \] con \(\hat{\xi}_i(\hat{s})\in (0,1)\).

En definitiva, podemos escribir: \[ \begin{align*} y_{i+1}-y_i & = h \sum_{k=0}^{m-1}C_k \nabla^{k} f(t_i,y(t_i))\\ & \quad +\frac{h^{m+1}}{m!}\int_{0}^{1}f^{(m)}(\hat{\xi}_i(\hat{s}),y(\hat{\xi}_i(\hat{s})))\hat{s} (\hat{s}+1)\cdots (\hat{s}+(m-1))\, d\hat{s}, \end{align*} \] con \(\displaystyle C_k = (-1)^k\int_0^1 \binom{-\hat{s}}{k}\,d\hat{s}\).

Hallemos los primeros \(C_k\): \[ \begin{align*} C_0 & = (-1)^0\int_0^1 \binom{-\hat{s}}{0}\,d\hat{s} =\int_0^1 1\, d\hat{s}=1,\\ C_1 & = (-1)^1\int_0^1 \binom{-\hat{s}}{1}\,d\hat{s} =-\int_0^1 -\hat{s}\, d\hat{s}=\left[\frac{\hat{s}^2}{2}\right]_0^1=\frac{1}{2},\\ C_2 & = (-1)^2\int_0^1 \binom{-\hat{s}}{2}\,d\hat{s} =\int_0^1 \frac{-\hat{s}(-\hat{s}-1)}{2}\, d\hat{s}=\frac{1}{2}\int_0^1 (\hat{s}^2+\hat{s})\,d\hat{s}\\ & =\frac{1}{2}\left[\frac{\hat{s}^3}{3}+\frac{\hat{s}^2}{2}\right]_0^1=\frac{1}{2}\left(\frac{1}{3}+\frac{1}{2}\right)=\frac{5}{12}, \end{align*} \]

\[ \begin{align*} C_3 & = (-1)^3\int_0^1 \binom{-\hat{s}}{3}\,d\hat{s} =-\int_0^1 \frac{-\hat{s}(-\hat{s}-1)(-\hat{s}-2)}{6}\, d\hat{s}\\ & =\frac{1}{6}\int_0^1 (\hat{s}^3+3\hat{s}^2+2\hat{s})\,d\hat{s} =\frac{1}{6}\left[\frac{\hat{s}^4}{4}+\hat{s}^3+\hat{s}^2\right]_0^1 \\ & =\frac{1}{6}\left(\frac{1}{4}+1+1\right)=\frac{3}{8}. \end{align*} \]

La tabla siguiente nos muestra los valores de los primeros \(C_k\):

A continuación, tratemos el término del error: \[ \frac{h^{m+1}}{m!}\int_{0}^{1}f^{(m)}(\hat{\xi}_i(\hat{s}),y(\hat{\xi}_i(\hat{s})))\hat{s} (\hat{s}+1)\cdots (\hat{s}+(m-1))\, d\hat{s}. \]

Teniendo en cuenta que la función \(\hat{s} (\hat{s}+1)\cdots (\hat{s}+(m-1))\) no cambia de signo en el intervalo \([0,1]\), usando el Teorema del valor medio para integrales podemos afirmar que existe un valor \(\hat{\xi}_i\in (0,1)\) tal que:

\[ \begin{align*} & \frac{h^{m+1}}{m!}\int_{0}^{1}f^{(m)}(\hat{\xi}_i(\hat{s}),y(\hat{\xi}_i(\hat{s})))\hat{s} (\hat{s}+1)\cdots (\hat{s}+(m-1))\, d\hat{s} \\ & =\frac{h^{m+1}f^{(m)}(\hat{\xi}_i,y(\hat{\xi}_i))}{m!}\int_{0}^{1}\hat{s} (\hat{s}+1)\cdots (\hat{s}+(m-1))\, d\hat{s}\\ & =h^{m+1} f^{(m)}(\hat{\xi}_i,y(\hat{\xi}_i))(-1)^m \int_{0}^{1} (-1)^m \binom{-\hat{s}}{m}\, d\hat{s}\\ & =h^{m+1} f^{(m)}(\hat{\xi}_i,y(\hat{\xi}_i))(-1)^m C_m. \end{align*} \]

En resumen, tenemos la expresión siguiente que nos relaciona \(y_{i+1}\) con \(y_i\): \[ y_{i+1}-y_i = h \sum_{k=0}^{m-1} C_k\nabla^{k} f(t_i,y(t_i))+h^{m+1} f^{(m)}(\hat{\xi}_i,y(\hat{\xi}_i))(-1)^m C_m. \]

La expresión anterior motiva la definición del método multipaso de orden \(m\) siguiente: \[ \hat{y}_{i+1} = \hat{y}_i+ h \sum_{k=0}^{m-1} C_k\nabla^{k} f(t_i,\hat{y}_i). \]

Dicho método se denomina método multipaso de Adams-Bashforth de orden \(m\) donde:

• Los coeficientes \(a_j\) valen \(a_0=\cdots =a_{m-2}=0\) y \(a_{m-1}=1\).
• Para hallar los coeficientes \(b_j\) hay que desarrollar la expresión \(\nabla^{k} f(t_i,y(t_i))\) y agrupar términos.

• Si escribimos la ecuación en diferencias que define el método multipaso de Adams-Bashforth de orden \(m\) de la forma: \[ \hat{y}_{i+1}=\hat{y}_i+h\phi_m(t_i,\hat{y}_i,\hat{y}_{i-1},\ldots,\hat{y}_{i+1-m}), \] con \(\displaystyle\phi_m(t_i,\hat{y}_i,\hat{y}_{i-1},\ldots,\hat{y}_{i+1-m}) = \sum_{k=0}^{m-1} C_k\nabla^{k} f(t_i,y(t_i))\), la relación que existe entre \(\phi_m\) y \(\phi_{m-1}\) es la siguiente: \[ \begin{align*} \phi_m(t_i,\hat{y}_i,\hat{y}_{i-1},\ldots,\hat{y}_{i+1-m})=& \phi_{m-1}(t_i,\hat{y}_i,\hat{y}_{i-1},\ldots,\hat{y}_{i+2-m})\\ & ^+C_{m-1}\nabla^{m-1} f(t_i,y(t_i)). \end{align*} \] Es decir, si tenemos calculado \(\phi_{m-1}(t_i,\hat{y}_i,\hat{y}_{i-1},\ldots,\hat{y}_{i+2-m})\), de cara a calcular \(\phi_m(t_i,\hat{y}_i,\hat{y}_{i-1},\ldots,\hat{y}_{i+1-m})\) basta calcular el término \(C_{m-1}\nabla^{m-1} f(t_i,y(t_i))\) y añadir dicho término a la expresión de \(\phi_{m-1}(t_i,\hat{y}_i,\hat{y}_{i-1},\ldots,\hat{y}_{i+2-m})\).

• El método multipaso anterior tiene orden de consistencia \(m\) ya que hemos probado que: \[ \begin{align*} & y_{i+1}-y_i = h \sum_{k=0}^{m-1} C_k\nabla^{k} f(t_i,y(t_i))+h^{m+1} f^{(m)}(\hat{\xi}_i,y(\hat{\xi}_i))(-1)^m C_m,\ \Rightarrow \\ &\frac{y_{i+1}-y_i}{h}-\sum_{k=0}^{m-1} C_k\nabla^{k} f(t_i,y(t_i)) = h^m f^{(m)}(\hat{\xi}_i,y(\hat{\xi}_i))(-1)^m C_m\\ & =O(h^m). \end{align*} \]

Calculemos la expresión de los métodos multipaso de Adams-Bashford para \(m=1,2,3,4\):

• \(m=1\): \[ \hat{y}_{i+1} = \hat{y}_i+ hf(t_i,\hat{y}_i), \] es el conocido método de Euler al ser un método de \(m=1\) paso.

• \(m=2\): \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i+ h \sum_{k=0}^{1} C_k\nabla^{k} f(t_i,\hat{y}_i) =\hat{y}_i +h\left(f(t_i,\hat{y}_i)+\frac{1}{2}\nabla f(t_i,\hat{y}_i)\right)\\ & = \hat{y}_i +h\left(f(t_i,\hat{y}_i)+\frac{1}{2}( f(t_i,\hat{y}_i)-f(t_{i-1},\hat{y}_{i-1}))\right)\\ & = \hat{y}_i + \frac{h}{2}(3 f(t_i,\hat{y}_i)-f(t_{i-1},\hat{y}_{i-1})). \end{align*} \]

• \(m=3\): \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i+ h \sum_{k=0}^{2} C_k\nabla^{k} f(t_i,\hat{y}_i) \\ & =\hat{y}_i + h\left(\frac{3}{2} f(t_i,\hat{y}_i)-\frac{1}{2}f(t_{i-1},\hat{y}_{i-1})+\frac{5}{12}\nabla^2 f(t_i,\hat{y}_i)\right)\\ & = \hat{y}_i + h\left(\frac{3}{2} f(t_i,\hat{y}_i)-\frac{1}{2}f(t_{i-1},\hat{y}_{i-1})+\right.\\ &\quad \left.+\frac{5}{12}(\nabla f(t_i,\hat{y}_i)-\nabla f(t_{i-1},\hat{y}_{i-1}))\right) \end{align*} \]

• \(m=3\) (continuación): \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i + h\left(\frac{3}{2} f(t_i,\hat{y}_i)-\frac{1}{2}f(t_{i-1},\hat{y}_{i-1})+\right.\\ &\quad \left.+\frac{5}{12}(f(t_i,\hat{y}_i)-f(t_{i-1},\hat{y}_{i-1})-(f(t_{i-1},\hat{y}_{i-1})-f(t_{i-2},\hat{y}_{i-2})))\right)\\ & = \hat{y}_i + h\left(\frac{23}{12}f(t_i,\hat{y}_i)-\frac{4}{3}f(t_{i-1},\hat{y}_{i-1})+\frac{5}{12}f(t_{i-2},\hat{y}_{i-2}\right)\\ & = \hat{y}_i+\frac{h}{12}(23 f(t_i,\hat{y}_i)-16 f(t_{i-1},\hat{y}_{i-1})+5 f(t_{i-2},\hat{y}_{i-2})). \end{align*} \]

• \(m=4\): \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i+ h \sum_{k=0}^{3} C_k\nabla^{k} f(t_i,\hat{y}_i) \\ & = \hat{y}_i+h\left(\frac{23}{12} f(t_i,\hat{y}_i)-\frac{16}{12} f(t_{i-1},\hat{y}_{i-1})+\frac{5}{12} f(t_{i-2},\hat{y}_{i-2})\right.\\ & \quad \left.+\frac{3}{8} \nabla^3 f(t_i,\hat{y}_i)\right)\\ & = \hat{y}_i+h\left(\frac{23}{12} f(t_i,\hat{y}_i)-\frac{16}{12} f(t_{i-1},\hat{y}_{i-1})+\frac{5}{12} f(t_{i-2},\hat{y}_{i-2})\right.\\ & \quad \left.+\frac{3}{8}(\nabla^2 f(t_i,\hat{y}_i)-\nabla^2 f(t_{i-1},\hat{y}_{i-1})\right) \end{align*} \]

• \(m=4\): \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i+h\Bigl(\frac{23}{12} f(t_i,\hat{y}_i)-\frac{16}{12} f(t_{i-1},\hat{y}_{i-1})+\frac{5}{12} f(t_{i-2},\hat{y}_{i-2})\\ & \quad +\frac{3}{8}(\nabla f(t_i,\hat{y}_i)-\nabla f(t_{i-1},\hat{y}_{i-1})\\ & \quad\quad -(\nabla f(t_{i-1},\hat{y}_{i-1})-\nabla f(t_{i-2},\hat{y}_{i-2}))\Bigr)\\ & = \hat{y}_i+h\Bigl(\frac{23}{12} f(t_i,\hat{y}_i)-\frac{16}{12} f(t_{i-1},\hat{y}_{i-1})+\frac{5}{12} f(t_{i-2},\hat{y}_{i-2})\\ & \quad +\frac{3}{8}(\nabla f(t_i,\hat{y}_i)-2\nabla f(t_{i-1},\hat{y}_{i-1})+\nabla f(t_{i-2},\hat{y}_{i-2})\Bigr) \end{align*} \]

• \(m=4\): \[ \begin{align*} \hat{y}_{i+1}& = \hat{y}_i+h\Bigl(\frac{23}{12} f(t_i,\hat{y}_i)-\frac{16}{12} f(t_{i-1},\hat{y}_{i-1})+\frac{5}{12} f(t_{i-2},\hat{y}_{i-2})\\ & \quad +\frac{3}{8}(f(t_i,\hat{y}_i)-f(t_{i-1},\hat{y}_{i-1})-2(f(t_{i-1},\hat{y}_{i-1})-f(t_{i-2},\hat{y}_{i-2}))\\ & \quad +f(t_{i-2},\hat{y}_{i-2})-f(t_{i-3},\hat{y}_{i-3})\Bigr) \\ & = \hat{y}_i +h\left(\frac{55}{24} f(t_i,\hat{y}_i)-\frac{59}{24}f(t_{i-1},\hat{y}_{i-1}+\frac{37}{24}f(t_{i-2},\hat{y}_{i-2})\right.\\ & \quad\left. -\frac{3}{8}f(t_{i-3},\hat{y}_{i-3})\right) \\ & = \hat{y}_i +\frac{h}{24}(55f(t_i,\hat{y}_i)-59 f(t_{i-1},\hat{y}_{i-1})+37 f(t_{i-2},\hat{y}_{i-2})\\ & \quad -9f(t_{i-3},\hat{y}_{i-3})) \end{align*} \]

• 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,2,3 (Calculamos los valores iniciales restantes \(\tilde{y}_1,\tilde{y}_2,\tilde{y}_3\) usando el método de RK-4. Supongamos que hemos definido la rutina RK4 de parámetros \(a\), \(b\), \(n\), \(y_0\), y \(f\), \(RK4(a,b,n,y_0,f)\))
  • Set \(t=t+h\)
  • Set \(\tilde{y}_i=RK4(a,t,i,y_0,f)\)
  • Set \(\mathbf{\hat{y}}_{i+1}=\tilde{y}_i\) (la componente \(i+1\)-ésima del vector de aproximaciones será la aproximación inicial \(i\)-ésima, \(i=1,2,3\))

• For i=4,...,n (Calculamos las demás aproximaciones)
  • Set \[ \begin{align*} \hat{y}_{i+1} & =\hat{y}_{i}+\frac{h}{24}(55\cdot f(t,\hat{y}_{i})-59\cdot f(t-h,\hat{y}_{i-1})\\ & \quad +37\cdot f(t-2\cdot h,\hat{y}_{i-2})-9f(t-3\cdot h,\hat{y}_{i-3})) \end{align*} \]
  • Set \(t=t+ h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos el vector de aproximaciones).

Recordemos que para deducir los métodos explícitos de Adams-Bashford, usábamos el polinomio de interpolación de la función \(f(t,y(t))\) en los nodos \(t_i,t_{i-1},\ldots,t_{i+1-m}\) para aproximar el valor de \(y(t_{i+1})=y_{i+1}\) en la integral de la expresión siguiente: \[ y_{i+1}=y_i+\int_{t_i}^{t_{i+1}}f(t,y(t))\,dt. \] Si añadimos el nodo \(t_{i+1}\) a los nodos anteriores \(t_i,t_{i-1},\ldots,t_{i+1-m}\), y repetimos todo el proceso realizado para deducir los métodos de Adams-Bashford, obtendremos los denominados métodos implícitos multipaso de Adams-Moulton de \(m\) pasos.

En primer lugar, podemos escribir la función \(f(t,y(t))\) como: \[ f(t,y(t))=P_{m+1}(t)+\frac{f^{(m+1)}(\xi_i(t),y(\xi_i(t)))}{(m+1)!}(t-t_{i+1})(t-t_{i})\cdots (t-t_{i+1-m}), \] donde \(P_{m+1}(t)\) es el polinomio interpolador de la función \(f(t,y(t))\) en los nodos equiespaciados \(t_{i+1},t_i,t_{i-1},\ldots,t_{i+1-m}\) que puede escribirse de la forma siguiente usando las diferencias hacia atrás: \[ P_{m+1}(t)=\sum_{k=0}^m (-1)^k \binom{s}{k}\nabla^k f(t_{i+1},y(t_{i+1})), \]

donde recordemos que \(s\) es la nueva variable que vale \(s=\frac{t_{i+1}-t}{h}\), \(\binom{s}{k}=\frac{s(s-1)\cdots (s-k+1)}{k!}\), y \[ \begin{align*} \nabla^0 f(t_{i+1},y(t_{i+1})) & =f(t_{i+1},y(t_{i+1})),\\ \nabla^k f(t_{i+1},y(t_{i+1})) & =\nabla^{k-1}f(t_{i+1},y(t_{i+1}))-\nabla^{k-1}f(t_{i},y(t_{i})). \end{align*} \]

Fijémonos que, como \(t\in [t_i,t_{i+1}]\), entonces \(s\in [0,1]\). Por tanto, haciendo el cambio de variable \(s=\frac{t_{i+1}-t}{h}\), \(dt=-hds\), tenemos que: \[ \begin{align*} y_{i+1}-y_i & = \int_{t_i}^{t_{i+1}}f(t,y(t))\, dt =h\int_0^{1}\sum_{k=0}^{m}(-1)^{k}\binom{s}{k}\nabla^{k} f(t_{i+1},y(t_{i+1}))\,ds\\ & \quad +\frac{h^{m+2}(-1)^{m+1}}{(m+1)!}\int_{0}^{1}f^{(m+1)}(\xi_i(s),y(\xi_i(s)))s (s-1)\cdots (s-m)\, ds, \end{align*} \] con \(\xi_i(s)\in (0,1)\).

En definitiva, podemos escribir: \[ \begin{align*} y_{i+1}-y_i & = h \sum_{k=0}^{m} \tilde{C}_k\nabla^{k} f(t_{i+1},y(t_{i+1}))\\ & \quad +\frac{h^{m+2}(-1)^{m+1}}{(m+1)!}\int_{0}^{1}f^{(m+1)}(\hat{\xi}_i(s),y(\hat{\xi}_i(s)))s (s-1)\cdots (s-m)\, ds, \end{align*} \] donde \(\displaystyle \tilde{C}_k = (-1)^k\int_0^1 \binom{s}{k}\,ds\).

Hallemos los primeros \(\tilde{C}_k\): \[ \begin{align*} \tilde{C}_0 & = (-1)^0\int_0^1 \binom{s}{0}\,ds =\int_0^1 1\, ds=1,\\ \tilde{C}_1 & = (-1)^1\int_0^1 \binom{s}{1}\,d\hat{s} =-\int_0^1 s\, ds=-\left[\frac{s^2}{2}\right]_0^1=-\frac{1}{2},\\ \tilde{C}_2 & = (-1)^2\int_0^1 \binom{s}{2}\,ds =\int_0^1 \frac{s(s-1)}{2}\, ds=\frac{1}{2}\int_0^1 (s^2-s)\,ds\\ & =\frac{1}{2}\left[\frac{s^3}{3}-\frac{s^2}{2}\right]_0^1=\frac{1}{2}\left(\frac{1}{3}-\frac{1}{2}\right)=-\frac{1}{12}, \end{align*} \]

\[ \begin{align*} \tilde{C}_3 & = (-1)^3\int_0^1 \binom{s}{3}\,d\hat{s} =-\int_0^1 \frac{s(s-1)(s-2)}{6}\, ds\\ & =-\frac{1}{6}\int_0^1 (s^3-3s^2+2s)\,ds =-\frac{1}{6}\left[\frac{s^4}{4}-s^3+s^2\right]_0^1 \\ & =-\frac{1}{6}\left(\frac{1}{4}-1+1\right)=-\frac{1}{24}. \end{align*} \]

La tabla siguiente nos muestra los valores de los primeros \(\tilde{C}_k\):

A continuación, tratemos el término del error: \[ \frac{h^{m+2}(-1)^{m+1}}{(m+1)!}\int_{0}^{1}f^{(m+1)}(\hat{\xi}_i(s),y(\hat{\xi}_i(s)))s (s-1)\cdots (s-m)\, ds. \]

Teniendo en cuenta que la función \(s (s-1)\cdots (s-m)\) no cambia de signo en el intervalo \([0,1]\), usando el Teorema del valor medio para integrales podemos afirmar que existe un valor \(\hat{\xi}_i\in (0,1)\) tal que:

\[ \begin{align*} & \frac{h^{m+2}(-1)^{m+1}}{(m+1)!}\int_{0}^{1}f^{(m+1)}(\hat{\xi}_i(s),y(\hat{\xi}_i(s)))s (s-1)\cdots (s-m)\, ds \\ & =h^{m+2}f^{(m+1)}(\xi_i,y(\xi_i))\int_{0}^{1}(-1)^{m+1}\binom{s}{m+1}\, ds\\ & =h^{m+2}f^{(m+1)}(\xi_i,y(\xi_i)) \tilde{C}_{m+1}. \end{align*} \]

En resumen, tenemos la expresión siguiente que nos relaciona \(y_{i+1}\) con \(y_i\): \[ y_{i+1}-y_i = h \sum_{k=0}^{m} \tilde{C}_k\nabla^{k} f(t_{i+1},y(t_{i+1}))+h^{m+2}f^{(m+1)}(\xi_i,y(\xi_i)) \tilde{C}_{m+1}. \]

La expresión anterior motiva la definición del método multipaso de orden \(m\) siguiente: \[ \hat{y}_{i+1} = \hat{y}_i+ h \sum_{k=0}^{m} \tilde{C}_k\nabla^{k} f(t_{i+1},\hat{y}_{i+1}). \]

Dicho método se denomina método multipaso de Adams-Moulton de orden \(m\) donde:

• Es un método implícito ya que \(\hat{y}_{i+1}\) aparece en los dos miembros de la ecuación en diferencias.
• Los coeficientes \(a_j\) valen \(a_0=\cdots =a_{m-2}=0\) y \(a_{m-1}=1\).
• Para hallar los coeficientes \(b_j\) hay que desarrollar la expresión \(\nabla^{k} f(t_{i+1},y(t_{i+1}))\) y agrupar términos.

• Si escribimos la ecuación en diferencias que define el método multipaso de Adams-Moulton de orden \(m\) de la forma: \[ \hat{y}_{i+1}=\hat{y}_i+h\tilde{\phi}_m(t_i,\hat{y}_{i+1},\hat{y}_i,\ldots,\hat{y}_{i+1-m}), \] con \(\displaystyle\tilde{\phi}_m(t_i,\hat{y}_{i+1},\hat{y}_i,\ldots,\hat{y}_{i+1-m}) = \sum_{k=0}^{m} \tilde{C}_k\nabla^{k} f(t_i,y(t_i))\), la relación que existe entre \(\phi_m\) y \(\phi_{m-1}\) es la siguiente: \[ \begin{align*} \tilde{\phi}_m(t_i,\hat{y}_{i+1},\hat{y}_i,\ldots,\hat{y}_{i+1-m})= & \tilde{\phi}_{m-1}(t_i,\hat{y}_{i+1},\hat{y}_i,\ldots,\hat{y}_{i+2-m})\\ & +C_m\nabla^{m} f(t_i,y(t_i)). \end{align*} \] Es decir, si tenemos calculado \(\tilde{\phi}_{m-1}(t_i,\hat{y}_{i+1},\hat{y}_i,\ldots,\hat{y}_{i+2-m})\), de cara a calcular \(\tilde{\phi}_m(t_i,\hat{y}_{i+1},\hat{y}_i,\ldots,\hat{y}_{i+1-m})\) basta calcular el término \(C_m\nabla^{m} f(t_i,y(t_i))\) y añadir dicho término a la expresión de \(\tilde{\phi}_{m-1}(t_i,\hat{y}_{i+1},\hat{y}_i,\ldots,\hat{y}_{i+2-m})\).

• El método multipaso anterior tiene orden de consistencia \(m+1\) ya que hemos probado que: \[ \begin{align*} & y_{i+1}-y_i = h \sum_{k=0}^{m} \tilde{C}_k\nabla^{k} f(t_{i+1},y(t_{i+1}))+h^{m+2}f^{(m+1)}(\xi_i,y(\xi_i)) \tilde{C}_{m+1},\ \Rightarrow \\ &\frac{y_{i+1}-y_i}{h}-\sum_{k=0}^{m} \tilde{C}_k\nabla^{k} f(t_{i+1},y(t_{i+1})) = h^{m+1}f^{(m+1)}(\xi_i,y(\xi_i)) \tilde{C}_{m+1}\\ & =O(h^{m+1}). \end{align*} \]

Calculemos la expresión de los métodos multipaso de Adams-Moulton para \(m=1,2,3\):

• \(m=1\): \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i+ h\sum_{k=0}^{1} \tilde{C}_k\nabla^{k} f(t_{i+1},\hat{y}_{i+1})\\ & = \hat{y}_i +h\left(f(t_{i+1},\hat{y}_{i+1})-\frac{1}{2}\nabla f(t_{i+1},\hat{y}_{i+1})\right)\\ & =\hat{y}_i +h\left(f(t_{i+1},\hat{y}_{i+1})-\frac{1}{2}(f(t_{i+1},\hat{y}_{i+1})-f(t_i,\hat{y}_i))\right)\\ & = \hat{y}_i+\frac{h}{2}(f(t_{i+1},\hat{y}_{i+1})+f(t_i,\hat{y}_i)), \end{align*} \] es el conocido método del trapecio al ser un método de \(m=1\) paso.

• \(m=2\): \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i+ h\sum_{k=0}^{2} \tilde{C}_k\nabla^{k} f(t_{i+1},\hat{y}_{i+1})\\ & =\hat{y}_i +h\left(\frac{1}{2}(f(t_{i+1},\hat{y}_{i+1})+f(t_{i},\hat{y}_{i})) -\frac{1}{12}\nabla^2 f(t_{i+1},\hat{y}_{i+1}) \right)\\ & =\hat{y}_i +h\left(\frac{1}{2}(f(t_{i+1},\hat{y}_{i+1})+f(t_{i},\hat{y}_{i})) \right.\\ & \quad \left. -\frac{1}{12}(\nabla f(t_{i+1},\hat{y}_{i+1})-\nabla f(t_i,\hat{y}_i))\right) \end{align*} \]

• \(m=2\) (continuación): \[ \begin{align*} \hat{y}_{i+1} & =\hat{y}_i +h\left(\frac{1}{2}(f(t_{i+1},\hat{y}_{i+1})+f(t_{i},\hat{y}_{i})) \right.\\ & \quad \left. -\frac{1}{12}(f(t_{i+1},\hat{y}_{i+1})-f(t_i,\hat{y}_i)-(f(t_i,\hat{y}_i)-f(t_{i-1},\hat{y}_{i-1})))\right) \\ & = \hat{y}_i + h \left(\frac{5}{12}f(t_{i+1},\hat{y}_{i+1})+\frac{2}{3}f(t_{i},\hat{y}_{i})-\frac{1}{12}f(t_{i-1},\hat{y}_{i-1})\right)\\ & = \hat{y}_i + \frac{h}{12}(5f(t_{i+1},\hat{y}_{i+1})+8f(t_{i},\hat{y}_{i})-f(t_{i-1},\hat{y}_{i-1})) \end{align*} \]

• \(m=3\): \[ \begin{align*} \hat{y}_{i+1} & = \hat{y}_i+ h \sum_{k=0}^{3} \tilde{C}_k\nabla^{k} f(t_{i+1},\hat{y}(t_{i+1})) \\ & =\hat{y}_i + h\left(\frac{5}{12}f(t_{i+1},\hat{y}_{i+1})+\frac{8}{12}f(t_{i},\hat{y}_{i})-\frac{1}{12}f(t_{i-1},\hat{y}_{i-1})\right.\\ & \quad \left.-\frac{1}{24}\nabla^3 f(t_{i+1},\hat{y}_{i+1}\right)\\ & =\hat{y}_i + h\left(\frac{5}{12}f(t_{i+1},\hat{y}_{i+1})+\frac{8}{12}f(t_{i},\hat{y}_{i})-\frac{1}{12}f(t_{i-1},\hat{y}_{i-1})\right.\\ & \quad \left.-\frac{1}{24} (\nabla^2 f(t_{i+1},\hat{y}_{i+1})-\nabla^2 f(t_{i},\hat{y}_{i})\right) \end{align*} \]

• \(m=3\) (continuación): \[ \begin{align*} \hat{y}_{i+1} & =\hat{y}_i + h\left(\frac{5}{12}f(t_{i+1},\hat{y}_{i+1})+\frac{8}{12}f(t_{i},\hat{y}_{i})-\frac{1}{12}f(t_{i-1},\hat{y}_{i-1})\right.\\ & \left.-\frac{1}{24} (\nabla f(t_{i+1},\hat{y}_{i+1})-\nabla f(t_{i},\hat{y}_{i})-(\nabla f(t_{i},\hat{y}_{i})-\nabla f(t_{i-1}-\hat{y}_{i-1}))\right)\\ & =\hat{y}_i + h\left(\frac{5}{12}f(t_{i+1},\hat{y}_{i+1})+\frac{8}{12}f(t_{i},\hat{y}_{i})-\frac{1}{12}f(t_{i-1},\hat{y}_{i-1})\right.\\ & \left.-\frac{1}{24} (\nabla f(t_{i+1},\hat{y}_{i+1})-2\nabla f(t_{i},\hat{y}_{i})+\nabla f(t_{i-1},\hat{y}_{i-1}))\right) \end{align*} \]

• \(m=3\) (continuación): \[ \begin{align*} \hat{y}_{i+1} & =\hat{y}_i + h\Bigl(\frac{5}{12}f(t_{i+1},\hat{y}_{i+1})+\frac{8}{12}f(t_{i},\hat{y}_{i})-\frac{1}{12}f(t_{i-1},\hat{y}_{i-1})\\ & -\frac{1}{24} (f(t_{i+1},\hat{y}_{i+1})-f(t_i,\hat{y}_i)-2(f(t_{i},\hat{y}_{i})-f(t_{i-1},\hat{y}_{i-1}))\\ &+ f(t_{i-1},\hat{y}_{i-1})-f(t_{i-2},\hat{y}_{i-2}))\Bigr)\\ & =\hat{y}_i +\frac{h}{24}(9f(t_{i+1},\hat{y}_{i+1})+19f(t_i,\hat{y}_i)-5f(t_{i-1},\hat{y}_{i-1})\\ & \quad +f(t_{i-2},\hat{y}_{i-2})). \end{align*} \]

• 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,2 (Calculamos los valores iniciales restantes \(\tilde{y}_1,\tilde{y}_2\) usando el método de RK-4. Supongamos que hemos definido la rutina RK4 de parámetros \(a\), \(b\), \(n\), \(y_0\), y \(f\), \(RK4(a,b,n,y_0,f)\))
  • Set \(t=t+h\)
  • Set \(\tilde{y}_i=RK4(a,t,i,y_0,f)\)
  • Set \(\mathbf{\hat{y}}_{i+1}=\tilde{y}_i\) (la componente \(i+1\)-ésima del vector de aproximaciones será la aproximación inicial \(i\)-ésima, \(i=1,2\))

• For i=3,...,n (Calculamos las demás aproximaciones)
  • Resolvemos la ecuación: \[ \begin{align*} \hat{y}_{i+1} & =\hat{y}_{i} +\frac{h}{24}(9f(t+h,\hat{y}_{i+1})+19f(t,\hat{y}_{i})-5f(t-h,\hat{y}_{i-1})\\ & \quad +f(t-2\cdot h,\hat{y}_{i-2})) \end{align*} \] usando un método numérico como Newton-Raphson, secante o regula-falsi.
  • Set \(t=t+ h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos el vector de aproximaciones).

Observación.

La ecuación a resolver en el paso \(i\)-ésimo es la siguiente: \[ \begin{align*} g(\hat{y}_{i+1}) & =\hat{y}_{i+1} -\hat{y}_{i} -\frac{h}{24}(9f(t+h,\hat{y}_{i+1})+19f(t,\hat{y}_{i})-5f(t-h,\hat{y}_{i-1})\\ & \quad +f(t-2\cdot h,\hat{y}_{i-2})) =0, \end{align*} \] donde recordemos que \(\hat{y}_{i+1}\) es la incógnita a calcular.

Si queremos usar el método de Newton-Raphson, tenemos que calcular la derivada de la función \(g(\hat{y}_{i+1})\): \[ g'(\hat{y}_{i+1})=1-\frac{3h}{8}\cdot \frac{\partial f}{\partial y}(t+h,\hat{y}_{i+1}). \]

Observación. (continuación)

Entonces, la sucesión de Newton-Raphson \((\hat{y}_{i+1})^{(n)}\) que converge a \(\hat{y}_{i+1}\) se tiene que calcular usando la recurrencia siguiente: \[ \hat{y}_{i+1}^{(n)}=\hat{y}_{i+1}^{(n-1)}-\frac{g(\hat{y}_{i+1}^{(n-1)})}{g'(\hat{y}_{i+1}^{(n-1)})}=\hat{y}_{i+1}^{(n-1)}-\frac{g(\hat{y}_{i+1}^{(n-1)})}{1-\frac{3h}{8}\cdot \frac{\partial f}{\partial y}(t+h,\hat{y}_{i+1}^{(n-1)})}, \] con \(\hat{y}_{i+1}^{(0)}=\hat{y}_{i}\).

En el ejemplo anterior, vimos que el método de Adams-Moulton de orden \(3\) con error de consistencia de orden \(4\), \(O(h^4)\) da mejores resultados que el método de Adams-Bashford de orden \(4\) con el mismo orden de consistencia.

En general, los métodos implícitos de Adams-Moulton tienen mejores resultados que los métodos explícitos de Adams-Bashsford.

Sin embargo, recordemos que los métodos implícitos son computacionalmente más costosos que los métodos explícitos al tener que hallar un cero de una función en cada paso de la iteración.

Por dicho motivo, los métodos implícitos son raramente usados ya que pesa más el coste computacional que la ganancia en los resultados obtenidos.

Una solución de compromiso es usar los denominados métodos predictor-corrector.

Éstos consisten en combinar un método explícito de Adams-Bashford al que llamamos predictor, \[ \hat{y}_{i+1}^{(p)}=\hat{y}_i^{(p)} + h\sum_{k=0}^{m-1}C_k\nabla^k f(t_i,\hat{y}_i^{(p)}), \] y un método implícito de Adams-Moulton, al que llamamos corrector, \[ \hat{y}_{i+1}^{(c)}=\hat{y}_i^{(c)} + h\sum_{k=0}^{m-1}\tilde{C}_k\nabla^k f(t_{i+1},\hat{y}_{i+1}^{(c)}), \] los dos del mismo orden de consistencia, \(O(h^{m})\) de la forma siguiente:

• Sea \(\hat{y}_i\) el valor de \(\hat{y}(t_i)\) en el paso \(t_i\).
• Para hallar \(\hat{y}_{i+1}\), el valor de \(\hat{y}(t_{i+1})\) en el paso \(t_{i+1}\),
  • en primer lugar calculamos \(\hat{y}_{i+1}^{(p)}\) usando la expresión del método predictor, \[ \hat{y}_{i+1}^{(p)}=\hat{y}_i+ h\sum_{k=0}^{m-1}C_k\nabla^k f(t_i,\hat{y}_i), \]
  • a continuación, sustituimos \(\hat{y}_{i+1}^{(c)}\) en el segundo miembro de la expresión del corrector por \(\hat{y}_{i+1}^{(p)}\) y hallamos \(\hat{y}_{i+1}\) usando su expresión: \[ \hat{y}_{i+1}=\hat{y}_i + h\sum_{k=0}^{m-1}\tilde{C}_k\nabla^k f(t_{i+1},\hat{y}_{i+1}^{(p)}). \]

Veamos cuál sería la expresión si usamos el método de Adams-Bashford de orden \(4\) y el método de Adams-Moulton de orden \(3\), los dos con orden de consistencia \(4\), \(O(h^4)\).

• Supongamos que tenemos calculado \(\hat{y}_i\).
• Para hallar \(\hat{y}_{i+1}\),
  • en primer lugar, hallamos \(\hat{y}_{i+1}^{(p)}\) usando el método de Adams-Bashford de orden \(4\): \[ \begin{align*} \hat{y}_{i+1}^{(p)} & =\hat{y}_i +\frac{h}{24}(55 f(t_i,\hat{y}_i)-59 f(t_{i-1},\hat{y}_{i-1})+37 f(t_{i-2},\hat{y}_{i-2})\\ & \quad -9 f(t_{i-3},\hat{y}_{i-3})), \end{align*} \]

• • por último, hallamos \(\hat{y}_{i+1}\) usando el método de Adams-Moulton de orden \(3\) sustituyendo \(\hat{y}_{i+1}^{(c)}\) por \(\hat{y}_{i+1}^{(p)}\): \[ \begin{align*} \hat{y}_{i+1} & =\hat{y}_i +\frac{h}{24}(9f(t_{i+1},\hat{y}_{i+1}^{(p)})+19 f(t_{i},\hat{y}_i)-5 f(t_{i-1},\hat{y}_{i-1})\\ & \quad ^+f(t_{i-2},\hat{y}_{i-2})). \end{align*} \]

Si iteramos el segundo paso, se pueden conseguir mejores aproximaciones \(\hat{y}_{i+1}\), pero dichas aproximaciones tienden al valor que obtendríamos si usamos el método de Adams-Moulton correspondiente.

• 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,2,3 (Calculamos los valores iniciales restantes \(\tilde{y}_1,\tilde{y}_2,\tilde{y}_3\) usando el método de RK-4. Supongamos que hemos definido la rutina RK4 de parámetros \(a\), \(b\), \(n\), \(y_0\), y \(f\), \(RK4(a,b,n,y_0,f)\))
  • Set \(t=t+h\)
  • Set \(\tilde{y}_i=RK4(a,t,i,y_0,f)\)
  • Set \(\mathbf{\hat{y}}_{i+1}=\tilde{y}_i\) (la componente \(i+1\)-ésima del vector de aproximaciones será la aproximación inicial \(i\)-ésima, \(i=1,2,3\))

• For i=4,...,n (Calculamos las demás aproximaciones)
  • Set \[ \begin{align*} \hat{y}_{i+1}^{(p)} & =\hat{y}_{i}+\frac{h}{24}(55\cdot f(t,\hat{y}_{i})-59\cdot f(t-h,\hat{y}_{i-1})\\ & \quad +37\cdot f(t-2\cdot h,\hat{y}_{i-2})-9f(t-3\cdot h,\hat{y}_{i-3})) \end{align*} \]
  • Set \[ \begin{align*} \hat{y}_{i+1} & =\hat{y}_{i} +\frac{h}{24}(9f(t+h,\hat{y}_{i+1}^{(p)})+19f(t,\hat{y}_{i})-5f(t-h,\hat{y}_{i-1})\\ & \quad +f(t-2\cdot h,\hat{y}_{i-2})) \end{align*} \]
  • Set \(t=t+ h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos el vector de aproximaciones).

Recordemos que un método multipaso para resolver el problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\) tenía como ecuación en diferencias:

\[ \begin{align*} \hat{y}_{i+1} & = a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m}+\\ &\quad h(b_m f(t_{i+1},\hat{y}_{i+1})+b_{m-1}f(t_i,\hat{y}_i)+\cdots +b_0 f(t_{i+1-m},\hat{y}_{i+1-m})), \end{align*} \] para \(i=m-1,m,\ldots,n-1\), donde se tienen que especificar los valores iniciales siguientes: \[ \hat{y}_0=y_0,\ \hat{y}_1=\tilde{y}_1,\ldots,\hat{y}_{m-1}=\tilde{y}_{m-1}. \]

La definición de convergencia para un método multipaso es la misma que para un método de un paso:

Definición de convergencia de un método multipaso.

Diremos que el método multipaso anterior es convergente si el error entre la solución aproximada y la solución exacta en los puntos del mallado \(t_i\) tiende a cero cuando \(h\) tiende a cero para todos los puntos de dicho mallado: \[ \lim_{h\to 0} |y(t_i)-\hat{y}_i| =0,\ i=0,\ldots,n. \]

Sin embargo, de cara a la consistencia, la definición es ligeramente diferente:

Definición de consistencia de un método multipaso.

Diremos que el método anterior es consistente si el error local de truncamiento \(\tau_{i+1}(h)\) tiende a cero si \(h\) tiende a cero: \[ \begin{align*} \lim_{h\to 0}\tau_{i+1}(h)= & \lim_{h\to 0}\frac{y(t_{i+1})-(a_{m-1}y_i +a_{m-2}y_{i-1}+\cdots +a_0 y_{i+1-m})}{h}\\ & -(b_m f(t_{i+1},\hat{y}_{i+1})+b_{m-1}f(t_i,\hat{y}_i)+\cdots +b_0 f(t_{i+1-m},\hat{y}_{i+1-m}))=0, \end{align*} \] \(i=0\ldots,n-1\), donde \(y_i=y(t_i)\) y además \(\displaystyle\lim_{h\to 0}|\tilde{y}_i-y(t_i)|=0,\ i=1,\ldots,m-1\).

El resultado siguiente nos dice cómo se comporta un método predictor-corrector con respecto a la consistencia y la convergencia:

Proposición. Consistencia y convergencia de un método predictor-corrector.

Consideremos un problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\) a la que hemos hallado una solución aproximada usando un método predictor de Adams-Bashford de orden \(m\) \[ \hat{y}_{i+1}=\hat{y}_i +h(b_{m-1}f(t_i,\hat{y}_i)+\cdots +b_0 f(t_{i+1-m},\hat{y}_{i+1-m})), \] con error local de truncamiento \(\tau_{i+1}^{(p)}(h)\) y un método corrector de Adams-Moulton de orden \(m-1\) \[ {\hat{y}}_{i+1}=\hat{y}_i +h(\tilde{b}_{m-1}f(t_{i+1},\hat{y}_{i+1})+\tilde{b}_{m-2}f(t_{i},\hat{y}_{i})+\cdots +\tilde{b}_0 f(t_{i+2-m},\hat{y}_{i+2-m})), \] con error local de truncamiento \(\tau_{i+1}^{(c)}(h)\).

Proposición. Consistencia y convergencia de un método predictor-corrector. (continuación)

Supongamos además que \(f(t,y)\) y \(\frac{\partial f}{\partial y}\) son continuas en \(D=\{(t,y),\ |\ a\leq t\leq b,\ -\infty <y<\infty\}\) y que \(\frac{\partial f}{\partial y}\) está acotada.

Entonces, el error local de truncamiento \(\tau_{i+1}(h)\) del método predictor-corrector es el siguiente: \[ \tau_{i+1}(h)=\tau_{i+1}^{(c)}(h)+\tau_{i+1}^{(p)}(h)\tilde{b}_{m-1}\frac{\partial f}{\partial y}(t_{i+1},\theta_{i+1}), \] donde \(\theta_{i+1}\in <0,h\tau_{i+1}^{(p)}(h)>\).

Proposición. Consistencia y convergencia de un método predictor-corrector. (continuación)

Además, existen constantes \(k_1\) y \(k_2\) tal que: \[ |\hat{y}_i-y_i|\leq \left(\max_{0\leq j\leq m-1}|\hat{y}_j-y_j|+k_1\tau(h)\right)\mathrm{e}^{k_2 (t_i-a)}, \] donde \(\tau(h)=\max_{m\leq j\leq n}|\tau_j(h)|.\)

Observación.

El resultado anterior nos dice que un método predictor-corrector con orden de consistencia \(m\) para los métodos predictor y corrector tiene orden de consistencia \(m\) y además el método es convergente, siempre y cuando el método para hallar las aproximaciones de las condiciones iniciales tenga también orden de consistencia \(m\).

Recordemos que un sistema de ecuaciones diferenciales de \(m\) ecuaciones con \(m\) incógnitas o \(m\) funciones reales de variable real es un sistema de ecuaciones del tipo: \[ \begin{align*} \frac{dy_1}{dt}=y_1'(t) & = f_1(t,y_1(t),\ldots,y_m(t)),\\ \frac{dy_2}{dt}=y_2'(t) & = f_2(t,y_1(t),\ldots,y_m(t)),\\ & \vdots\\ \frac{dy_m}{dt}=y_m'(t) & = f_m(t,y_1(t),\ldots,y_m(t)), \end{align*} \] donde ahora las “incógnitas” son las funciones \(y_1(t),y_2(t),\ldots,y_m(t)\), para \(t\in [a,b]\) y condiciones iniciales: \[ y_1(a)=y_{01},\ y_2(a)=y_{02},\ldots,y_m(a)=y_{0m}. \]

El objetivo es hallar unas soluciones aproximadas \(\hat{y}_1(t),\hat{y}_2(t),\ldots,\hat{y}_m(t)\).

Los valores en los que se hallarán dichas soluciones aproximadas serán en valores discretos de \(t\).

Es decir, fijado un valor \(h>0\), consideramos el mallado siguiente \(t_i=a+i\cdot h\), donde \(i=0,\ldots,n\) y por tanto, la relación entre \(h\) y \(n\) tiene que ser \(h=\frac{b-a}{n}\) para que el último valor del mallado sea \(t_n=b\), el valor final de \(t\) en el intervalo de integración.

Entonces las soluciones aproximadas se hallarán en los valores de mallado \(t_i\): \[ \hat{y}_1(t_i),\hat{y}_2(t_i),\ldots,\hat{y}_m(t_i),\ i=0,\ldots,m, \] con \[ \hat{y}_1(a)=y_{01},\ \hat{y}_2(a)=y_{02},\ldots,\hat{y}_m(a)=y_{0m}. \]

Para que un sistema de ecuaciones diferenciales tenga solución, necesitamos que las funciones \(f_i\) verifiquen la condición de Lipschitz respecto las \(m\) últimas variables:

Definición de condición de Lipschitz.

Sea \(f(t,y_1,\ldots,y_m)\) una función de \(m+1\) variables definida en el dominio: \[ D=\{(t,u_1,\ldots,u_m),\ |\ a\leq t\leq b,\ -\infty <u_i <\infty,\ i=1,2,\ldots,m\}, \] satisface la condición de Lipschitz en \(D\) respecto las variables \(u_1,\ldots,u_m\) si existe una constante \(L\), denominada constante de Lipschitz de \(f\) tal que para todo \((t,u_1,\ldots,u_m), (t,u_1',\ldots,u_m')\in D\), \[ |f(t,u_1,\ldots,u_m)-f(t,u_1',\ldots,u_m')|\leq L\sum_{i=1}^m |u_i-u_i'|. \]

Verificar la condición de Lipschitz de una función \(m+1\) variables es bastante tedioso a partir de la definición.

El siguiente resultado nos simplifica las cosas para funciones de clase \({\cal C}^1\):

Proposición. Condición suficiente de condición de Lipschitz.

Sea \(f(t,y_1,\ldots,y_m)\) una función de \(m+1\) variables de clase \({\cal C}^1\) definida en el dominio: \[ D=\{(t,u_1,\ldots,u_m),\ |\ a\leq t\leq b,\ -\infty <u_i <\infty,\ i=1,2,\ldots,m\}. \] Supongamos que existe una constante \(L\) tal que para cualquier valor \((t,u_1,\ldots,u_m)\in D\), se cumple \(\left|\frac{\partial f}{\partial u_i}(t,u_1,\ldots,u_m)\right|\leq L\), \(i=1,\ldots,m\).

Entonces \(f\) verifica la condición de Lipschitz en el dominio \(D\) con constante de Lipschitz \(L\).

El resultado siguiente nos dice cuando un sistema de ecuaciones diferenciales tiene solución única:

Teorema de existencia y unicidad de soluciones.

Sea \(D\) el dominio siguiente: \[ D=\{(t,u_1,\ldots,u_m),\ |\ a\leq t\leq b,\ -\infty <u_i <\infty,\ i=1,2,\ldots,m\}. \] Supongamos que las funciones \(f_i(t,u_1,\ldots,u_m)\) son continuas y verifican la condición de Lipschitz con constante de Lipschitz \(L_i\).

Teorema de existencia y unicidad de soluciones (continuación).

Entonces el sistema de ecuaciones diferenciales \[ \begin{align*} \frac{dy_1}{dt}=y_1'(t) & = f_1(t,y_1(t),\ldots,y_m(t)),\\ \frac{dy_2}{dt}=y_2'(t) & = f_2(t,y_1(t),\ldots,y_m(t)),\\ & \vdots\\ \frac{dy_m}{dt}=y_m'(t) & = f_m(t,y_1(t),\ldots,y_m(t)), \end{align*} \] con condiciones iniciales: \[ y_1(a)=y_{01},\ y_2(a)=y_{02},\ldots,y_m(a)=y_{0m}, \] tiene solución única \(y_1(t),\ldots,y_m(t)\), \(a\leq t\leq b\).

Observación.

Considerando \(\displaystyle L=\max_{i=1,\ldots,m}L_i\), podemos afirmar que todas las funciones \(f_i(t,u_1,\ldots,u_m)\) verifican la condición de Lipschitz con la misma constante de Lipschitz \(L\).

Todos los métodos numéricos para resolver el problema de valores iniciales correspondiente a la ecuación diferencial \(y'=f(t,y)\), \(a\leq t\leq b\), \(y(a)=y_0\) se pueden adaptar al problema de valores iniciales correspondiente al sistema de ecuaciones diferenciales: \[ \begin{align*} \frac{dy_1}{dt}=y_1'(t) & = f_1(t,y_1(t),\ldots,y_m(t)),\\ \frac{dy_2}{dt}=y_2'(t) & = f_2(t,y_1(t),\ldots,y_m(t)),\\ & \vdots\\ \frac{dy_m}{dt}=y_m'(t) & = f_m(t,y_1(t),\ldots,y_m(t)), \end{align*} \] con condiciones iniciales: \[ y_1(a)=y_{01},\ y_2(a)=y_{02},\ldots,y_m(a)=y_{0m}. \]

Para no complicar la notación, supondremos que los métodos numéricos son explícitos.

Supongamos que la ecuación en diferencias del método numérico para resolver la ecuación diferencial es de la forma: \[ \hat{y}_{i+1}=\hat{y}_i +h\varphi_f(t_i,\hat{y}_i),\ \hat{y}_0=y_0,\ i=0,\ldots,n-1, \] donde hemos escrito \(\varphi_f\) indicando la dependencia de la función \(\varphi_f\) respecto la función \(f\).

Entonces, dicho método aplicado al sistema de ecuaciones diferenciales constará del siguiente sistema en diferencias:

\[ \begin{align*} \hat{y}_{i+1}^{(1)} & =\hat{y}_i^{(1)} +h \varphi_{f_1}(t_i,\hat{y}_i^{(1)},\hat{y}_i^{(2)},\ldots,\hat{y}_i^{(m)}),\\ \hat{y}_{i+1}^{(2)} & =\hat{y}_i^{(2)} +h \varphi_{f_2}(t_i,\hat{y}_i^{(1)},\hat{y}_i^{(2)},\ldots,\hat{y}_i^{(m)}),\\ & \vdots\\ \hat{y}_{i+1}^{(m)} & =\hat{y}_i^{(m)} +h \varphi_{f_m}(t_i,\hat{y}_i^{(1)},\hat{y}_i^{(2)},\ldots,\hat{y}_i^{(m)}),\\ \end{align*} \] con \(\hat{y}_0^{(1)}=y_{01},\ \hat{y}_0^{(2)}=y_{02},\ldots,\hat{y}_0^{(m)}=y_{0m}\), con \(i=0,\ldots,n-1\).

Recordemos la ecuación en diferencias del método de Euler: \[ \hat{y}_{i+1} = \hat{y}_i +h f(t_i,\hat{y}_i),\ i=0,\ldots,n-1,\ \hat{y}_0=y_0. \] En este caso, tenemos que \(\varphi_f(t,\hat{y})=f(t,\hat{y})\). Por tanto, el sistema en diferencias será en este caso: \[ \begin{align*} \hat{y}_{i+1}^{(1)} & =\hat{y}_i^{(1)} +h f_1(t_i,\hat{y}_i^{(1)},\hat{y}_i^{(2)},\ldots,\hat{y}_i^{(m)}),\\ \hat{y}_{i+1}^{(2)} & =\hat{y}_i^{(2)} +h f_2(t_i,\hat{y}_i^{(1)},\hat{y}_i^{(2)},\ldots,\hat{y}_i^{(m)}),\\ & \vdots\\ \hat{y}_{i+1}^{(m)} & =\hat{y}_i^{(m)} +h f_m(t_i,\hat{y}_i^{(1)},\hat{y}_i^{(2)},\ldots,\hat{y}_i^{(m)}),\\ \end{align*} \] con \(\hat{y}_0^{(1)}=y_{01},\ \hat{y}_0^{(2)}=y_{02},\ldots,\hat{y}_0^{(m)}=y_{0m}\), con \(i=0,\ldots,n-1\).

• INPUT valores inicial y final \(a\), \(b\), entero \(n\), condición inicial \(\mathbf{y}_0=(y_{01},\ldots,y_{0m})\) y funciones \(f_1(t,y_1,\ldots,y_m),\ldots,f_m(t,y_1,\ldots,y_m)\)
• Set \(h=\frac{b-a}{n}\).
• Set \(t=a\).
• Set \(\mathbf{\hat{y}}=\mathbf{0}\). (Definimos inicialmente la matriz de aproximaciones de dimensiones \(m\times (n+1)\) como \(0\). Vamos a guardar en la columna \(i\)-ésima de la matriz \(\mathbf{\hat{y}}\) la aproximación \(\mathbf{\hat{y}}(t+(i-1)h)\))
• Set \(\mathbf{\hat{y}}(,1)=\mathbf{y}_0\). (Definimos la primera columna de la matriz de aproximaciones como las condiciones iniciales)

• For i=1,...,n
  • For j=1,...,m
    • Set \(\mathbf{\hat{y}}(j,i+1)=\mathbf{\hat{y}}(j,i)+h\cdot f_j(t,\mathbf{\hat{y}}(1,i),\ldots,\mathbf{\hat{y}}(m,i))\). (calculamos la componente \(j\)-ésima de \(\mathbf{\hat{y}}(t+ih)\))
  • Set \(t=t+h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos la matriz de aproximaciones).

Recordemos que la idea de los métodos de Taylor era desarrollar el valor \(\hat{y}(t_{i+1})=\hat{y}_{i+1}\) alrededor de \(t=t_i\) usando el desarrollo de Taylor pero llegando hasta orden \(k\) donde si \(k=1\), tenemos el método de Euler: \[ \hat{y}_{i+1}= \hat{y}_i+h \hat{y}'(t_i)+\frac{h^2}{2}\hat{y}''(t_i)+\cdots +\frac{h^k}{k!}\hat{y}^{(k)}(t_i). \] El método de Taylor de orden 2 consiste en desarrollar el desarrollo anterior para \(k=2\): \[ \hat{y}_{i+1}= \hat{y}_i+h \hat{y}'(t_i)+\frac{h^2}{2}\hat{y}''(t_i). \] para resolver el problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\).

En el caso de sistemas, tenemos que aplicar el desarrollo anterior para cada componente \(\hat{y}_i^{(j)}\), \(j=1,\ldots,m\):

\[ \hat{y}_{i+1}^{(j)}= \hat{y}_i^{(j)}+h \hat{y}'(t_i)^{(j)}+\frac{h^2}{2}\hat{y}''(t_i)^{(j)}. \]

El valor de \(\hat{y}'(t_i)^{(j)}\) será el siguiente usando que \(\hat{y}'(t_i)\approx f_j(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)})\): \[ \hat{y}'(t_i)^{(j)}= f_j(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)}). \]

Para calcular \(\hat{y}''(t_i)^{(j)}\) usaremos la regla de la cadena: \[ \begin{align*} \hat{y}''(t_i)^{(j)}\approx & f_j'(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)})\\ = & \frac{\partial f_j}{\partial t}(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)})+\sum_{l=1}^m\frac{\partial f_j}{\partial y_l}(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)})\cdot y_l'(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)})\\ = & \frac{\partial f_j}{\partial t}(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)})+\sum_{l=1}^m\frac{\partial f_j}{\partial y_l}(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)})\cdot f_l(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)}). \end{align*} \]

Usando la notación siguiente introducida anteriormente, \[ \begin{align*} f_{i,j}= & f_j(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)}),\\ f_{t,i,j}^{(1)}=&\frac{\partial f_j}{\partial t}(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)}),\\ f_{y,i,l,j}^{(1)}=&\frac{\partial f_j}{\partial y_l}(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)}), \end{align*} \]

el desarrollo anterior será: \[ \begin{align*} \hat{y}_{i+1}^{(j)}=& \hat{y}_i^{(j)}+h f_{i,j}+\frac{h^2}{2}\left(f_{t,i,j}^{(1)}+\sum_{l=1}^m f_{y,i,l,j}^{(1)}\cdot f_{i,l}\right)\\ = & \hat{y}_i^{(j)}+h \left(f_{i,j}+\frac{h}{2}\left(f_{t,i,j}^{(1)}+\sum_{l=1}^m f_{y,i,l,j}^{(1)}\cdot f_{i,l}\right)\right), \end{align*} \] \(j=1,\ldots,m\), \(i=0,\ldots,n-1\).

Entonces, dicho método aplicado al sistema de ecuaciones diferenciales constará del siguiente sistema en diferencias:

\[ \begin{align*} \hat{y}_{i+1}^{(1)} & =\hat{y}_i^{(1)} +h \left(f_{i,1}+\frac{h}{2}\left(f_{t,i,1}^{(1)}+\sum_{l=1}^m f_{y,i,l,1}^{(1)}\cdot f_{i,l}\right)\right),\\ \hat{y}_{i+1}^{(2)} & =\hat{y}_i^{(2)} +h \left(f_{i,2}+\frac{h}{2}\left(f_{t,i,2}^{(1)}+\sum_{l=1}^m f_{y,i,l,2}^{(1)}\cdot f_{i,l}\right)\right),\\ & \vdots\\ \hat{y}_{i+1}^{(m)} & =\hat{y}_i^{(m)} +h \left(f_{i,m}+\frac{h}{2}\left(f_{t,i,m}^{(1)}+\sum_{l=1}^m f_{y,i,l,m}^{(1)}\cdot f_{i,l}\right)\right), \end{align*} \] con \(\hat{y}_0^{(1)}=y_{01},\ \hat{y}_0^{(2)}=y_{02},\ldots,\hat{y}_0^{(m)}=y_{0m}\), con \(i=0,\ldots,n-1\).

• INPUT valores inicial y final \(a\), \(b\), entero \(n\), condición inicial \(\mathbf{y}_0=(y_{01},\ldots,y_{0m})\) y funciones \(f_j(t,y_1,\ldots,y_m)\), \(j=1,\ldots,m\), \(f_{t,j}^{(1)}(t,y_1,\ldots,y_m)=\frac{\partial f_j}{\partial t}(t,y_1,\ldots,y_m)\), \(j=1,\ldots,m\) y \(f_{y,l,j}^{(1)}(t,y_1,\ldots,y_m)=\frac{\partial f_j}{\partial y_l}(t,y_1,\ldots,y_m)\), \(l,j=1,\ldots,m\).
• Set \(h=\frac{b-a}{n}\).
• Set \(t=a\).
• Set \(\mathbf{\hat{y}}=\mathbf{0}\). (Definimos inicialmente la matriz de aproximaciones de dimensiones \(m\times (n+1)\) como \(0\). Vamos a guardar en la columna \(i\)-ésima de la matriz \(\mathbf{\hat{y}}\) la aproximación \(\mathbf{\hat{y}}(t+(i-1)h)\))
• Set \(\mathbf{\hat{y}}(,1)=\mathbf{y}_0\). (Definimos la primera columna de la matriz de aproximaciones como las condiciones iniciales)

• For i=1,...,n
  • For j=1,...,m
    • Set \[ \begin{align*} & \mathbf{\hat{y}}(j,i+1)= \mathbf{\hat{y}}(j,i)+h\cdot (f_{j}(t,\mathbf{\hat{y}}(1,i),\ldots,\mathbf{\hat{y}}(m,i))\\ & +(h/2)\cdot (f_{t,j}^{(1)}(t,\mathbf{\hat{y}}(1,i),\ldots,\mathbf{\hat{y}}(m,i))\\ & +\sum_{l=1}^m f_{y,l}^{(1)}(t,\mathbf{\hat{y}}(1,i),\ldots,\mathbf{\hat{y}}(m,i)) \cdot f_l(t,\mathbf{\hat{y}}(1,i),\ldots,\mathbf{\hat{y}}(m,i)))). \end{align*} \] (calculamos la componente \(j\)-ésima de \(\mathbf{\hat{y}}(t+ih)\))
  • Set \(t=t+h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos la matriz de aproximaciones).

Recordemos que un método de Runge-Kutta explícito general es un método de un paso donde la expresión de la ecuación en diferencias era la 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_{p=1}^J c_p k_p,\\ k_p & = f\left(t_i+h a_p,\hat{y}_i+h\sum_{l=1}^{p-1} b_{pl} k_l\right),\ p=1,\ldots,J. \end{align*} \] donde \(\displaystyle\sum_{k=1}^{p-1} b_{pl}=a_p\).

En el caso de sistemas, la generalización de la expresión anterior para la componente \(\hat{y}_i^{(j)}\), \(j=1,\ldots,m\) sería: \(\hat{y}_{i+1}^{(j)}=\hat{y}_i^{(j)}+h\varphi_j(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)}),\) donde \[ \begin{align*} \varphi_j(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)}) & =\sum_{p=1}^J c_p k_p^{(j)},\\ k_p^{(j)} & = f_j\left(t_i+h a_p,\hat{y}_i^{(1)}+h\sum_{l=1}^{p-1} b_{pl} k_l^{(1)},\ldots, \hat{y}_i^{(m)}+h\sum_{l=1}^{p-1} b_{pl} k_l^{(m)}\right), \end{align*} \] \(p=1,\ldots,J,\ j=1,\ldots,m.\)

Desarrollemos la expresión anterior para el caso de Runge-Kutta 4.

Recordemos la expresión de la ecuación en diferencias finitas para resolver el problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\): \[ \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*} \]

Las ecuaciones en diferencias finitas para sistemas serían las siguientes: \[ \begin{align*} \hat{y}_{i+1}^{(j)} & = \hat{y}_i^{(j)} +\frac{h}{6}(k_1^{(j)}+2k_2^{(j)}+2k_3^{(j)}+k_4^{(j)}),\\ k_1^{(j)} & = f_j(t_i,\hat{y}_i^{(1)},\ldots,\hat{y}_i^{(m)}),\\ k_2^{(j)} & = f_j\left(t_i+\frac{h}{2},\hat{y}_i^{(1)}+\frac{h k_1^{(1)}}{2},\ldots,\hat{y}_i^{(m)}+\frac{h k_1^{(m)}}{2}\right),\\ k_3^{(j)} & = f_j\left(t_i+\frac{h}{2},\hat{y}_i^{(1)}+\frac{h k_2^{(1)}}{2},\ldots,\hat{y}_i^{(m)}+\frac{h k_2^{(m)}}{2}\right),\\ k_4^{(j)} & = f_j(t_i+h,\hat{y}_i^{(1)} + h k_3^{(1)},\ldots,\hat{y}_i^{(m)} + h k_3^{(m)}),\ \hat{y}_0^{(1)}=y_{01},\ldots,\hat{y}_{0m}=y_{0m}, \end{align*} \] \(i=0,\ldots,n-1\), \(j=1,\ldots,m\).

• INPUT valores inicial y final \(a\), \(b\), entero \(n\), condición inicial \(\mathbf{y}_0=(y_{01},\ldots,y_{0m})\) y funciones \(f_j(t,y_1,\ldots,y_m)\), \(j=1,\ldots,m\).
• Set \(h=\frac{b-a}{n}\).
• Set \(t=a\).
• Set \(\mathbf{\hat{y}}=\mathbf{0}\). (Definimos inicialmente la matriz de aproximaciones de dimensiones \(m\times (n+1)\) como \(0\). Vamos a guardar en la columna \(i\)-ésima de la matriz \(\mathbf{\hat{y}}\) la aproximación \(\mathbf{\hat{y}}(t+(i-1)h)\))
• Set \(\mathbf{\hat{y}}(,1)=\mathbf{y}_0\). (Definimos la primera columna de la matriz de aproximaciones como las condiciones iniciales)

• For i=1,...,n
• Set \(\mathbf{k}=\mathbf{0}\) (En la matriz \(\mathbf{k}\) de dimensiones \(m\times 4\) guardaremos \(k_1^{(j)},k_2^{(j)},k_3^{(j)},k_4^{(j)},\ j=1,\ldots,m\), la columna \(l\)-ésima de la matriz \(\mathbf{k}\) será \(k_l^{(1)},\ldots,k_l^{(m)}\), \(l=1,2,3,4\))
  • For j=1,...,m
    • Set \(\mathbf{k}(j,1)=f_j(t,\mathbf{\hat{y}}(1,i),\ldots,\mathbf{\hat{y}}(m,i))\)
  • For j=1,...,m
    • Set \(\mathbf{k}(j,2)=f_j\left(t+\frac{h}{2},\mathbf{\hat{y}}(1,i)+\frac{h\cdot k(1,1)}{2},\ldots,\mathbf{\hat{y}}(m,i)+\frac{h\cdot k(m,1)}{2}\right)\)

• • For j=1,...,m
    • Set \(\mathbf{k}(j,3)=f_j\left(t+\frac{h}{2},\mathbf{\hat{y}}(1,i)+\frac{h\cdot k(1,2)}{2},\ldots,\mathbf{\hat{y}}(m,i)+\frac{h\cdot k(m,2)}{2}\right)\)
  • For j=1,...,m
    • Set \(\mathbf{k}(j,4)=f_j\left(t+h,\mathbf{\hat{y}}(1,i)+h\cdot k(1,3),\ldots,\mathbf{\hat{y}}(m,i)+h\cdot k(m,3)\right)\)
  • For j=1,...,m
    • Set \(\mathbf{\hat{y}}(i+1,j)=\mathbf{\hat{y}}(i,j)+\frac{h}{6}(\mathbf{k}(j,1)+2\cdot\mathbf{k}(j,2)+2\cdot\mathbf{k}(j,3)+\mathbf{k}(j,4))\). (calculamos la componente \(j\)-ésima de \(\mathbf{\hat{y}}(t+ih)\))
  • Set \(t=t+h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos la matriz de aproximaciones).

Un problema de valores iniciales asociado a una ecuación diferencial de orden \(m\) consiste en hallar una solución aproximada \(\hat{y}(t)\) de la ecuación diferencial siguiente: \[ y^{(m)}(t)=f(t,y, y',y'',\ldots,y^{(m-1)}),\ a\leq t\leq b, \] con condiciones iniciales \(y(a)=y_0,\ y'(a)=y'_0,\ y''(a)=y''_0,\ldots, y^{(m-1)}(a)=y^{(m-1)}_0\), donde \(y^{(j)}\) indica la derivada \(j\)-ésima de la función \(y(t)\) y los valores \(y_0,y'_0,\ldots,y^{(m-1)}_0\) se suponen conocidos.

Veamos que resolver un problema de valores iniciales asociado a una ecuación diferencial de orden \(m\) es equivalente a resolver un problema de valores iniciales asociado a un sistema de ecuaciones diferenciales de \(m\) ecuaciones.

Para ver dicha equivalencia, dado una ecuación diferencial de orden \(m\): \[ y^{(m)}(t)=f(t,y, y',y'',\ldots,y^{(m-1)}),\ a\leq t\leq b, \] consideramos las funciones siguientes asociadas a dicha ecuación diferencial: \[ y_1(t)=y(t),\ y_2(t)=y'(t),\ldots,y_m(t)=y^{(m-1)}(t),\ a\leq t\leq b. \]

Usando las funciones anteriores, tendremos que el sistema diferencial que verifican dichas funciones es el siguiente: \[ \begin{align*} \frac{dy_1}{dt}= & y_1'(t) = f_1(t,y_1(t),\ldots,y_m(t))=y_2(t),\\ \frac{dy_2}{dt}= & y_2'(t) = f_2(t,y_1(t),\ldots,y_m(t))=y_3(t),\\ & \vdots\\ \frac{dy_{m-1}}{dt}= & y_{m-1}'(t) = f_{m-1}(t,y_1(t),\ldots,y_m(t))=y_{m}(t),\\ \frac{dy_m}{dt}= & y_m'(t) = f_m(t,y_1(t),\ldots,y_m(t))=y^{(m)}(t)=f(t,y_1(t),y_2(t),\ldots,y_{m-1}(t)), \end{align*} \] con \(a\leq t\leq b\) con condiciones iniciales \(y_1(a)=y_0,\ y_2(a)=y'(a)=y'_0,\ y_3(a)=y''(a)=y''_0,\ldots,\) \(y_{m}(a)=y^{(m-1)}(a)=y^{(m-1)}_0\).

Es decir, las funciones \(f_i(t,y_1(t),\ldots,y_m(t))\), \(i=1,\ldots,m\) son las siguientes en el caso de ecuaciones diferenciales de orden \(m\): \[ \begin{align*} f_1(t,y_1(t),\ldots,y_m(t))= &y_2(t),\\ f_2(t,y_1(t),\ldots,y_m(t))= & y_3(t),\\ & \vdots\\ f_{m-1}(t,y_1(t),\ldots,y_m(t))= & y_{m}(t),\\ f_m(t,y_1(t),\ldots,y_m(t))= & f(t,y_1(t),y_2(t),\ldots,y_{m-1}(t)), \end{align*} \]

Observación.

Recordemos que para que un sistema de ecuaciones diferenciales tenga solución única, las funciones \(f_i(t,y_1(t),\ldots,y_m(t))\) tenían que verifican la condición de Lipschitz con constante de Lipschitz \(L_i\).

En el caso de ecuaciones diferenciales de orden \(m\), después de considerar el cambio al sistema diferencial correspondiente de orden \(m\), las funciones \(f_i(t,y_1(t),\ldots,y_m(t))=y_{i+1}(t)\), \(i=1,\ldots,m-1\) verifican la condición de Lipschitz con constante de Lipschitz \(L_i=1\) ya que: \[ |f_i(t,u_1,\ldots,u_m)-f_i(t,u_1',\ldots,u_m')|=|u_{i+1}-u_{i+1}'|\leq\sum_{j=1}^m |u_j-u_j'|. \]

En resumen, para que una ecuación diferencial de orden \(m\) tenga solución única, es suficiente que la función \(f_m(t,y_1(t),\ldots,y_{m}(t))=f(t,y(t),y'(t),\ldots,y^{(m-1)}(t))\) verifique la condición de Lipschitz con constante de Lipschitz \(L_m\).