Las ecuaciones diferenciales ordinarias son una herramienta fundamental en la resolución de problemas en ingeniería, machine learning y en ciencia en general.

Una ecuación diferencial es una ecuación del tipo: \[ \frac{dy}{dt}=y'(t)=f(t,y(t)),\ \mathrm{para\ }a\leq t\leq b, \] donde \(t\) es la variable independiente (en algunos libros la indican por \(x\)) que suponemos que está en un intervalo determinado \([a,b]\) e \(y(t)\) es la función que tenemos que hallar que verifique la ecuación anterior.

Por tanto, la “incógnita” de una ecuación diferencial es una función real de variable real indicada por \(y(t)\) de la que tenemos que hallar su valor en un intervalo determinado \([a,b]\).

En caso en que hubiese una solución de la ecuación anterior, ésta no es única. Para “imponer” la unicidad, tenemos que dar una condición inicial del tipo \(y(a)=y_0\).

Entonces resolver la ecuación diferencial anterior, con la condición inicial \(y(a)=y_0\) se denomina resolver un problema de valor inicial para ecuaciones diferenciales.

En este capítulo, trataremos también el caso de tener un sistema de ecuaciones diferenciales. 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}. \]

Una ecuación diferencial del tipo \(\frac{dy}{dt}=y'(t)=f(t,y)\) se denomina ecuación diferencial de primer orden ya que sólo aparece la derivada primera.

Entonces una ecuación diferencial de orden \(n\) es una ecuación del tipo: \[ y^{(n)}=f(t,y(t),y'(t),y''(t),\ldots,y^{(n-1)}(t)), \] para \(t\in [a,b]\) en un cierto intervalo \([a,b]\) con condiciones iniciales: \[ y(a)=y_{0},\ y'(a)=y_1,\ldots y^{(n-1)}(a)=y_{n-1}. \] Veremos cómo transformar una ecuación diferencial de orden \(n\) en un sistema de ecuaciones diferenciales con \(n\) ecuaciones.

Entonces resolver una ecuación diferencial de orden \(n\) es equivalente a hallar la solución de un sistema de \(n\) ecuaciones diferenciales.

Las ecuaciones diferenciales aparecen en problemas en ingeniería y ciencia cuando intentamos modelizar un fenómeno que queremos estudiar.

Una vez establecida la ecuación diferencial existen dos perspectivas para hallar una solución aproximada del fenómeno:

• Simplificar la ecuación diferencial por otra de la que seamos capaces de hallar una solución analítica exacta y usar dicha solución como solución aproximada de nuestra ecuación diferencial o
• Usar métodos específicos de aproximación de la ecuación diferencial original y será la perspectiva que desarrollaremos en este capítulo.

Antes de pasar a estudiar cómo hallar una solución aproximada de una ecuación diferencial, recordemos los conceptos necesarios para saber cuándo dicha ecuación tiene solución.

Dada una ecuación diferencial de la forma \(y'=f(t,y)\), con \(t\in [a,b]\) en un cierto intervalo \([a,b]\) con condición inicial \(y(a)=y_0\), para estudiar si dicha ecuación tiene solución, necesitamos que la función de dos variables \(f\) verifique la condición de Lipschitz:

Definición de condición de Lipschitz. Dada una función \(f: D\longrightarrow \mathbb{R}\) de dos variables \(f(t,y)\) definida en un dominio \(D\subset\mathbb{R}^2\), diremos que verifica la condición de Lipschitz respecto la segunda variable \(y\) si existe una constante \(L>0\) tal que para cualquier \((t,y_1)\in D\), \((t,y_2)\in D\), se cumple: \[ |f(t,y_1)-f(t,y_2)|\leq L |y_1-y_2|. \]

La constante \(L\) se denomina constante de Lipschitz de la función \(f\).

Otra de las condiciones que necesitamos es que el dominio de definición de la función \(f\) sea convexo:

Definición de convexidad. Sea \(D\subset \mathbb{R}^2\) un subconjunto del plano \(\mathbb{R}^2\). Diremos que el conjunto \(D\) es convexo si para cualquier par de puntos \((t_1,y_1),(t_2,y_2)\in D\) del conjunto \(D\), el segmento que une los dos puntos anteriores también está en \(D\), es decir, para cualquier \(\lambda\in [0,1]\), \[ \lambda (t_1,y_1)+(1-\lambda) (t_2,y_2)=(\lambda t_1+(1-\lambda)t_2,\lambda y_1+(1-\lambda)y_2)\in D. \]

Comprobar que un función \(f\) verifica la condición de Lipschitz es bastante tedioso si usamos la definición. Sin embargo, el teorema siguiente nos puede facilitar las cosas:

Teorema. Condición suficiente de Lipschitz. Dada una función \(f: D\longrightarrow \mathbb{R}\) de dos variables \(f(t,y)\) de clase \({\cal C}^1\) definida en un dominio \(D\subset\mathbb{R}^2\) convexo tal que para todo \((t,y)\in D\) se cumple que existe una constante \(L\) tal que: \[ \left|\frac{\partial f}{\partial y}(t,y)\right|\leq L. \] Entonces la función \(f\) verifica la condición de Lipschitz respecto la segunda variable \(y\) con constante de Lipschitz \(L\).

Resumamos lo visto hasta el momento:

• Dada una ecuación diferencial \(y'=f(t,y)\) con \(t\in [a,b]\) con condición inicial \(y(a)=y_0\), sabemos que tiene solución única si la función \(f\) verifica la condición de Lipschitz con constante de Lipschitz \(L\).

• Además, si el dominio de definición \(D\) es convexo, para verificar que la función \(f\) verifica la condición de Lipschitz, basta ver que la parcial con respecto \(y\), \(\frac{\partial f}{\partial y}\), está acotada en el dominio \(D\) por una constante \(L\) y dicha constante sería la constante de Lipschitz.

Ahora bien, como nuestro objetivo es hallar una solución aproximada de la ecuación diferencial, no nos conformamos en que exista dicha solución, sino que necesitamos que dicha solución sea estable, es decir, que si tenemos pequeños errores de redondeo en la condición inicial, \(y_0\) o en el cálculo de la función \(f(t,y(t))\), la solución aproximada \(\hat{y}(t)\) esté controlada o tenga errores que estén controlados de alguna manera.

Definamos formalmente lo dicho anteriormente:

Definición de problema bien planteado. Diremos que el problema de valor inicial \(\frac{dy}{dt}=y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0,\) es un problema bien planteado si:

• existe una única solución \(y(t)\) y
• existen constantes \(\epsilon_0\) y \(C>0\) tal que para cualquier \(\epsilon\in (0,\epsilon_0)\) y para cualquier función \(\delta (t)\) continua, con \(|\delta(t)|<\epsilon\), para todo \(t\in [a,b]\) y si el valor \(\delta_0\) verifica \(|\delta_0| <\epsilon\), el problema de valor inicial perturbado: \[ \frac{dz}{dt}=z'=f(t,z)+\delta(t),\ a\leq t\leq b,\ z(a)=y_0+\delta_0, \] tiene solución única \(z(t)\) que verifica: \[ |z(t)-y(t)|\leq C\epsilon,\ \mathrm{para\ todo\ }t\in [a,b]. \]

Observación.

La definición anterior nos dice que si perturbamos la condición inicial con una constante \(\delta_0\) y la ecuación diferencial con una función \(\delta(t)\), ambas acotadas por un valor \(\epsilon\) que se supone pequeño, entonces la solución del problema perturbado \(z(t)\) dista de la solución del problema no perturbado \(y(t)\) como máximo en \(C\epsilon\), es decir, está controlada y no puede tener un error tan grande como se quiera.

Ver que un problema está bien planteado a partir de la definición es un trabajo arduo muy complicado. Sin embargo, el teorema siguiente nos simplifica las cosas:

Teorema. Condición suficiente de problema bien planteado.

Supongamos que el dominio \(D\) es de la forma: \[ D=\{(t,y),\ a\leq t\leq b,\ -\infty<y<\infty\}. \] Sea una función \(f: D\longrightarrow \mathbb{R}\) de dos variables \(f(t,y)\) continua que satisface la condición de Lipschitz respecto la segunda variable \(y\) de constante \(L\). Entonces el problema de condición inicial: \[\frac{dy}{dt}=y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0,\] es un problema bien planteado.

Observación.

El teorema anterior nos dice que si \(D\) es un rectángulo “infinito” y por tanto convexo y \(f\) una función continua que verifica la condición de Lipschitz respecto la segunda variable \(y\), el problema automáticamente está bien planteado y serán el tipo de problemas que analizaremos en este capítulo.

Vamos a empezar a ver cómo calcular una solución aproximada del problema de valor inicial \(y'=f(t,y)\) con \(t\in [a,b]\) y \(y(a)=y_0\).

Notaremos con \(\hat{y}\) a la solución aproximada del problema anterior.

No calcularemos \(\hat{y}\) en todos los valores \(t\in [a,b]\) sino en un conjunto discreto de \(t\)’s. Más concretamente,

• en primer lugar fijaremos un paso \(h>0\) y seguidamente
• iremos calculando \(\hat{y}\) en los valores \(t_i=a+ih\), para \(i=0,\ldots,n\).

Para asegurarnos que “llegamos” al valor final \(t=b\), la relación entre \(h\) y \(n\) debe ser la siguiente: \(h=\frac{b-a}{n}\) ya que de esta forma, \(t_n=a+nh=a+b-a=b\) y calcularemos \(\hat{y}(t_n)=\hat{y}(b)\).

El valor de \(\hat{y}(a)\) será el valor inicial y no tiene error \(\hat{y}(a)=y(a)=y_0\).

Los demás valores \(\hat{y}(t_i)=\hat{y}(a+ih)\) se irán calculando en función de \(h\), los valores anteriores \(\hat{y}(t_j)\), y de los valores \(f(t_j,\hat{y}(t_j))\), \(j=0,\ldots,i-1\).

Iremos desarrollando métodos de tal forma que la solución aproximada \(\hat{y}\) y la solución \(y\) sean los más “próximas” posible en el mallado \(t_i=a+ih\), \(i=0,\ldots,n\).

La definición de proximidad en estos momentos es un concepto vago que iremos desarrollando y puliendo en las siguientes subsecciones.

Una vez establecidas las nociones básicas necesarias para calcular la solución aproximada del problema de valor inicial \(\frac{dy}{dt}=y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0,\) pasemos a estudiar los métodos más conocidos para hallar dicha solución aproximada \(\hat{y}\).

Los métodos más sencillos son los denominados métodos de un paso.

Una vez fijado el paso \(h\) y el mallado \(t_i=a+ih\), \(i=0,\ldots,n\), un método de un paso tiene la forma: \[ \hat{y}(t_{i+1})=\hat{y}(t_i)+h\varphi(t_i,\hat{y}(t_i)),\ i=0,\ldots,n-1, \] Entonces los métodos de un paso vienen determinados por la función \(\varphi\).

Se llaman métodos de un paso ya que los valores de la solución aproximada en \(t_{i+1}\) se obtienen a partir de los valores de la solución aproximada en el valor anterior \(t_i\) del mallado con \(\hat{y}(t_0)=\hat{y}(a)=y_0\).

De cara a simplificar la notación anterior, llamaremos \(\hat{y}_i\) al valor de la solución aproximada \(\hat{y}\) en el punto del mallado \(t_i\): \(\hat{y}_i=\hat{y}(t_i)\).

De esta forma, el método de un paso se escribe de la forma siguiente: \[ \hat{y}_{i+1}=\hat{y}_i + h\varphi(t_i,\hat{y}_i),\ i=0,\ldots,n-1,\ \hat{y}_0=y_0. \] Así, los demás valores de la solución aproximada se calcularán de la forma siguiente: \[ \begin{align*} \hat{y}_1 & = \hat{y}_0 +h\varphi(t_0,\hat{y}_0)=y_0 +h\varphi(a,y_0),\\ \hat{y}_2 & = \hat{y}_1 +h\varphi(t_1,\hat{y}_1),\\ & \vdots\\ \hat{y}_n & =\hat{y}(b) = \hat{y}_{n-1} +h\varphi(t_{n-1},\hat{y}_{n-1}). \end{align*} \]

Analicemos el método de un paso más sencillo: el conocido método de Euler.

Dicho método aparece cuando desarrollamos el valor de la solución aproximada en el punto del mallado \(t_{i+1}=t_i+h\), \(\hat{y}_{i+1}=\hat{y}(t_i+h)\), usando la siguiente expresión del desarrollo de Taylor de la función \(\hat{y}(t)\): \[ \hat{y}_{i+1}=\hat{y}(t_i+h)\approx \hat{y}(t_i)+h\cdot \hat{y}'(t_i)=y_i +h\cdot \hat{y}'(t_i). \] Ahora bien, como \(y'=f(t,y)\) y \(\hat{y}\) es una solución aproximada, podemos escribir que \(\hat{y}'(t_i)=f(t_i,\hat{y}(t_i))=f(t_i,\hat{y}_i)\).

Entonces, la expresión anterior será: \[ \hat{y}_{i+1} = \hat{y}_i +h f(t_i,\hat{y}_i),\ i=0,\ldots,n-1,\ \hat{y}_0=y_0. \] Por tanto, el método de Euler es un método de un paso donde la función \(\varphi\) vale: \(\varphi(t,\hat{y})=f(t,\hat{y})\).

• 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\cdot f(t,\hat{y}_{i})\). (calculamos \(\hat{y}(t+i h)\))
  • Set \(t=t+h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos el vector de aproximaciones).

El teorema siguiente nos da la cota del error para el método de Euler:

Teorema. Cota del error para el método de Euler.

Sea \(f(t,y)\) una función de dos variables continua y que cumple la condición de Lipschitz con respecto a la segunda variable con constante de Lipschitz \(L\) en la región: \[ D=\{(t,y),\ a\leq t\leq b,\ -\infty<y<\infty\}. \] Sea \(y(t)\) la solución única del problema de valores iniciales: \[ \frac{dy}{dt}=y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0. \] Supongamos que la función \(y(t)\) es de clase \({\cal C}^2\) y que su derivada segunda está acotada en el intervalo \([a,b]\), es decir, existe una constante \(M\) tal que \(|y''(t)|\leq M\), para todo \(t\in [a,b]\).

Teorema. Cota del error para el método de Euler (continuación).

Sea \(n>0\), un entero positivo, \(h=\frac{b-a}{n}\) y \(t_i=a+ih\), \(i=0,\ldots,n\) el mallado correspondiente al paso \(h\).

Sea \(\hat{y}(t_i)=\hat{y}_i\) la solución aproximada por el método de Euler en el mallado anterior.

Entonces: \[ |y(t_i)-\hat{y}(t_i)|=|y(t_i)-\hat{y}_i|\leq \frac{hM}{2L}\left(\mathrm{e}^{L(t_i-a)}-1\right)=\frac{hM}{2L}\left(\mathrm{e}^{iLh}-1\right), \] para \(i=0,1,\ldots,n\).

Para demostrar el teorema, necesitamos dos lemas previos:

Lema 1.

Sea \(x\geq -1\) un valor real mayor o igual que \(-1\). Entonces para cualquier valor \(m>0\), se cumple que \[ 0\leq (1+x)^m\leq \mathrm{e}^{mx}. \]

Lema 2.

Sean \(s>0\), \(t>0\) dos valores reales positivos. Sean \((a_n)_{n=0}^N\) una sucesión finita de números reales verificando que \[ a_0\geq -\frac{s}{t},\quad a_{n+1}\leq (1+s)a_n+t,\ n=0,1,\ldots,N-1. \] Entonces: \[ a_{n+1}\leq \mathrm{e}^{(n+1)s}\left(a_0+\frac{t}{s}\right)-\frac{t}{s}. \]

Observación.

En el ejemplo anterior hemos usado que conocemos el valor de la solución \(y(t)\) para hallar la cota de la derivada segunda \(M\).

En general, no hace falta conocer la expresión de la solución exacta \(y(t)\) de cara a calcular una expresión de \(y''(t)\) ya que ésta se puede calcular de la forma siguiente usando la regla de la cadena: \[ \begin{align*} y'(t) & =f(t,y(t)),\ \Rightarrow \\ y''(t) & =\frac{\partial f}{\partial t}(t,y(t))+\frac{\partial f}{\partial y}(t,y(t))\cdot y'(t) =\frac{\partial f}{\partial t}(t,y(t))+\frac{\partial f}{\partial y}(t,y(t))\cdot f(t,y(t)), \end{align*} \] pudiéndose obtener una expresión para \(y''(t)\).

Dado un método numérico de un paso \(\hat{y}_{i+1} =\hat{y}_{i}+h\varphi(t_{i},\hat{y}_{i})\) con \(\hat{y}_0=y_0\), necesitamos una manera de estimar el error que cometemos respecto la solución exacta \(y(t)\) del problema de valor inicial \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\).

Dado un valor del mallado \(t_i=a+ih\), un primer intento sería estimar el error cometido \(|y(t_i)-\hat{y}_i|\) pero pensemos que en la mayoría de los casos desconocemos la función que nos da la solución exacta \(y(t_i)\) haciendo impracticable calcular una cota de dicho error.

Para intentar solventar el problema anterior, se introduce el denominado error local de truncamiento o LTE en sus siglas en inglés: Local Truncation Error.

La solución aproximada \(\hat{y}_i\) es un vector de \(n\) componentes real que soluciona la ecuación en diferencias dada por el método numérico en cuestión: \[ \hat{y}_{i+1}=\hat{y}_i+h \varphi(t_i,\hat{y}_i),\ i=0,\ldots,n-1,\ \hat{y}_0=y_0. \]

Entonces el error local de truncamiento, \(\tau_{i+1}(h)\), es el error cometido por la solución exacta \(y(t_i)\) al intentar verificar la ecuación en diferencias anterior: \[ \begin{align*} \tau_{i+1}(h) & =\frac{y(t_{i+1})-(y(t_i)+h\varphi(t_i,y(t_i)))}{h}\\ & =\frac{y(t_{i+1})-y(t_i)}{h}-\varphi(t_i,y(t_i)),\ i=0,\ldots,n-1. \end{align*} \]

Podemos interpretar el error local de truncamiento como la diferencia entre el cociente incremental de la solución exacta \(y(t_i)\) usando diferencias hacia adelante (ver la primera parte del curso de métodos numéricos) y la función \(\varphi(t_i,y(t_i))\) que nos caracteriza el método numérico en cuestión.

La idea es: cuanto mejor aproxima la función \(\varphi\) del método numérico el cociente incremental de la aproximación de la derivada \(y'(t)\) usando la solución exacta, mejor será el método numérico en cuestión.

Dado un método numérico, si demostramos que \(\tau_{i+1}(h)=O(h^m)\), diremos que el método numérico tiene orden \(m\).

Cuánto más grande sea el valor \(m\), mejor será el método numérico. De esta forma, el error local de truncamiento sirve para comparar lo “buenos” que son los métodos numéricos considerados.

Veamos que el método de Euler tiene un error local de truncamiento de orden \(1\), es decir, es un método muy limitado.

Recordemos que la ecuación en diferencias que verificaba el vector de aproximaciones \(\hat{y}_i\), \(i=0,\ldots,n-1\) era: \[ \hat{y}_{i+1}=\hat{y}_i+h f(t_i,\hat{y}_i),\ \hat{y}_0=y_0. \]

Entonces el error local de truncamiento será: \[ \begin{align*} \tau_{i+1}(h) & =\frac{y(t_{i+1})-y(t_i)}{h}-f(t_i,y(t_i))=\frac{y(t_i+h)-y(t_i)}{h}-f(t_i,y(t_i)) \\ & =\frac{y(t_i)+hy'(t_i)+\frac{h^2}{2}y''(\xi_i)-y(t_i)}{h}-f(t_i,y(t_i))\\ & = y'(t_i)+\frac{h}{2}y''(\xi_i)-f(t_i,y(t_i)), \end{align*} \] con \(\xi_i\in (t_i,t_{i+1})\). Usando que \(y(t)\) es la solución exacta de la ecuación diferencial, podemos afirmar que \(y'(t_i)=f(t_i,y(t_i))\). Por tanto, \[ \tau_{i+1}(h) = \frac{h}{2}y''(\xi_i) =O(h). \] En conclusión, el método de Euler tiene error local de truncamiento de orden \(1\).

El error local de truncamiento motiva la definición siguiente:

Definición de consistencia de un método.

Consideramos el método numérico siguiente: \[ \hat{y}_{i+1}=\hat{y}_i+h\varphi(t_i,\hat{y}_i),\ i=0,\ldots,n-1,\ \hat{y}_0=y_0, \] para resolver el problema de valor inicial \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\).

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: \(\displaystyle\lim_{h\to 0}\tau_{i+1}(h)=0.\)

Además, diremos que el método es consistente de orden \(k\) si el error local de truncamiento \(\tau_{i+1}(h)\) verifica \(\displaystyle\tau_{i+1}(h)=O(h^k)\).

Observación.

Si un método es consistente de orden \(k\) con \(k\geq 1\), entonces el método es consistente.

La defininción de convergencia de un método es la siguiente:

Definición de convergencia de un método.

Consideramos el método numérico siguiente: \[ \hat{y}_{i+1}=\hat{y}_i+h\varphi(t_i,\hat{y}_i),\ i=0,\ldots,n-1,\ \hat{y}_0=y_0, \] para resolver el problema de valor inicial \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\).

Diremos que el método 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. \]

Además, diremos que un método es convergente de orden \(k\) si para cualquier punto del mallado \(t_i\) se cumple: \(\displaystyle |y(t_i)-\hat{y}_i|=O(h^k)\).

Observación.

La última condición de convergencia puede expresarse de la forma siguiente de forma equivalente: \[ \lim_{h\to 0} \max_{i=0,\ldots,n} |y(t_i)-\hat{y}_i| =0. \]

Vamos a generalizar el método de Euler. Recordemos que introducíamos el método de Euler desarrollando por Taylor el valor \(\hat{y}(t_{i+1})=\hat{y}_{i+1}\) alrededor de \(t=t_i\) hasta llegar a orden \(1\).

La idea será pues 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). \]

Los valores \(y^{(j)}(t_i)\), usando que \(\hat{y}\) es una aproximación de la ecuación diferencial \(y'=f(t,y)\), se calculan de la forma siguiente: \[ \hat{y}'(t_i)=f(t_i,\hat{y}_i),\ \hat{y}''(t_i)=f'(t_i,\hat{y}_i),\ldots,\hat{y}^{(j)}(t_i)=f^{(j-1)}(t_i,\hat{y}_i), \] donde \(f^{(j-1)}(t_i,\hat{y}_i)\) es la derivada \(j-1\)-ésima de la función real de variable real \(f(t,y(t))\) con respecto la variable \(t\) en \(t=t_i\).

Calculemos los primeros valores usando la regla de la cadena: \[ \begin{align*} \hat{y}'(t_i) & =f(t_i,\hat{y}_i),\\ \hat{y}''(t_i) & =f'(t_i,\hat{y}_i) =\frac{\partial f}{\partial t}(t_i,\hat{y}_i)+\frac{\partial f}{\partial y}(t_i,\hat{y}_i)\hat{y}'(t_i)\\ & =\frac{\partial f}{\partial t}(t_i,\hat{y}_i)+\frac{\partial f}{\partial y}(t_i,\hat{y}_i)f(t_i,\hat{y}_i), \end{align*} \]

\[ \begin{align*} \hat{y}'''(t_i) & =f''(t_i,\hat{y}_i) =\frac{\partial^2 f}{\partial t^2}(t_i,\hat{y}_i)+\frac{\partial^2 f}{\partial t\partial y}(t_i,\hat{y}_i)f(t_i,\hat{y}_i)+\frac{\partial^2 f}{\partial y\partial t}(t_i,\hat{y}_i)f(t_i,\hat{y}_i)\\ & \quad+\frac{\partial^2 f}{\partial y^2}(t_i,\hat{y}_i)f^2(t_i,\hat{y}_i)+ \frac{\partial f}{\partial y}(t_i,\hat{y}_i)\left(\frac{\partial f}{\partial t}(t_i,\hat{y}_i)+\frac{\partial f}{\partial y}(t_i,\hat{y}_i)f(t_i,\hat{y}_i)\right)\\ & = \frac{\partial^2 f}{\partial t^2}(t_i,\hat{y}_i)+2\frac{\partial^2 f}{\partial t\partial y}(t_i,\hat{y}_i)f(t_i,\hat{y}_i)+\frac{\partial^2 f}{\partial y^2}(t_i,\hat{y}_i)f^2(t_i,\hat{y}_i)\\ & \quad +\frac{\partial f}{\partial y}(t_i,\hat{y}_i)\frac{\partial f}{\partial t}(t_i,\hat{y}_i)+\left(\frac{\partial f}{\partial y}(t_i,\hat{y}_i)\right)^2 f(t_i,\hat{y}_i). \end{align*} \]

Como vemos que la notación hace tedioso escribir las derivadas \(f^{(j)}(t_i,\hat{y}_i)\) vamos a simplificarla de la forma siguiente: \[ \begin{align*} f(t_i,\hat{y}_i) & =f_i,\\ \frac{\partial^j f}{\partial t^j}(t_i,\hat{y}_i) & =f_{t,i}^{(j)},\\ \frac{\partial^j f}{\partial y^j}(t_i,\hat{y}_i)& =f_{y,i}^{(j)},\\ \frac{\partial^{j+l} f}{\partial t^j\partial y^l}(t_i,\hat{y}_i) & =f_{t,y,i}^{(j,l)} \end{align*} \]

Entonces las derivadas anteriores se escribirán de la forma siguiente: \[ \begin{align*} \hat{y}'(t_i) & =f_i,\\ \hat{y}''(t_i) & =f'(t_i,\hat{y}_i) =f_{t,i}^{(1)}+f_{y,i}^{(1)}f_i,\\ \hat{y}'''(t_i) & =f''(t_i,\hat{y}_i) = f_{t,i}^{(2)}+2 f_{t,y,i}^{(1,1)}f_i+f_{y,i}^{(2)}f_i^2+f_{y,i}^{(1)}f_{t,i}^{(1)}+\left(f_{y,i}^{(1)}\right)^2 f_i,\\ \hat{y}^{(iv)}(t_i) & = f'''(t_i,\hat{y}_i) = f_{y,i}^{(3)} f_i^3+\left(4 f_{y,i}^{(1)} f_{y,i}^{(2)} +3 f_{t,y,i}^{(1,2)}\right)f_i^2\\ &\quad +\left((f_{y,i}^{(1)})^3+5 f_{t,y,i}^{(1,1)} f_{y,i}^{(1)}+3\left(f_{y,i}^{(2)} f_{t,i}^{(1)} + f_{t,y,i}^{(2,1)}\right)\right) f_i\\ & \quad +f_{t,i}^{(1)}\left((f_{y,i}^{(1)})^2+3 f_{t,y,i}^{(1,1)}\right)+f_{y,i}^{(1)} f_{t,i}^{(2)}+f_{t,i}^{(3)} \\ & \vdots \end{align*} \]

El método de Taylor de orden 2 resuelve la ecuación en diferencias siguiente: \[ \hat{y}_{i+1}=\hat{y}_i + h\left(f_i + \frac{h}{2} (f_{t,i}^{(1)}+f_{y,i}^{(1)} f_i)\right),\ \hat{y}_0=y_0, \] donde recordemos que: \[ f_i=f(t_i,\hat{y}_i),\quad f_{t,i}^{(1)}=\frac{\partial f}{\partial t}(t_i,\hat{y}_i),\quad f_{y,i}^{(1)}=\frac{\partial f}{\partial y}(t_i,\hat{y}_i) \]

• INPUT valores inicial y final \(a\), \(b\), entero \(n\), condición inicial \(y_0\) y funciones \(f(t,y)\), \(f_t(t,y)=\frac{\partial f}{\partial t}(t,y)\) y \(f_y(t,y)=\frac{\partial f}{\partial y}(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\cdot (f(t,\hat{y}_{i})+\frac{h}{2}\cdot (f_t(t,\hat{y}_{i})+f_y(t,\hat{y}_{i})\cdot f(t,\hat{y}_{i})))\). (calculamos \(\hat{y}(t+i h)\))
  • Set \(t=t+h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos el vector de aproximaciones).

El siguiente resultado nos dice que el método de Taylor orden \(k\) es consistente de orden \(k\):

Error local de truncamiento en los métodos de Taylor.

Consideremos el siguiente problema de valores iniciales: \[ y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0, \] donde suponemos que \(f\) es de Lipschitz en la segunda variable y el dominio de definición de \(f\) es de la forma \[ D=\{(t,y)\in\mathbb{R}^2, \ a\leq t\leq b,\ -\infty < y <\infty\}. \] Suponemos que la función \(f\) es de clase \({\cal C}^k\). Entonces el método de Taylor de orden \(k\) para hallar una solución aproximada al problema de valores iniciales anterior es consistente de orden \(k\).

Observación.

Como \(f\) es de clase \({\cal C}^k\) tendremos que la solución \(y(t)\) que verifica el problema de valores iniciales será de clase \({\cal C}^{k+1}\) ya que \(y^{(k+1)}(t)=f^{(k)}(t,y(t))\).

Hasta ahora, hemos supuesto que la función \(\varphi\) del método numérico de un paso: \[ \hat{y}_{i+1} =\hat{y}_i+h\varphi(t_i,\hat{y}_i), \] sólo depende de \(t_i\) y de \(\hat{y}_i\).

A este tipo de métodos se les denomina explícitos ya que podemos calcular el valor de la solución aproximada \(\hat{y}_{i+1}\) en el valor de mallado \(t_{i+1}\) directamente a partir del valor de la solución aproximada \(\hat{y}_i\) en el valor del mallado \(t_i\).

Ahora bien, si suponemos que la función \(\varphi\) depende de \(t_i,\hat{y}_i\) y \(\hat{y}_{i+1}\), tendremos el método de un paso siguiente: \[ \hat{y}_{i+1} =\hat{y}_i+h\varphi(t_i,\hat{y}_i,\hat{y}_{i+1}). \] En este caso, no podemos calcular directamente \(\hat{y}_{i+1}\) en función de \(\hat{y}_i\) y de \(t_i\) ya que la “incógnita” \(\hat{y}_{i+1}\) aparece en los dos miembros de la definición del método.

Para hallar dicho valor \(\hat{y}_{i+1}\), tendremos que hallar un cero de la ecuación siguiente: \[ \hat{y}_{i+1} -\hat{y}_i- h\varphi(t_i,\hat{y}_i,\hat{y}_{i+1})=0. \] A este tipo de métodos se les denomina implícitos.

Como podemos observar, en general, los métodos implícitos son computacionalmente más costosos que los explícitos.

Sin embargo, con expresiones de la misma complejidad que los métodos explícitos, pueden alcanzar un orden de consistencia mayor.

Veamos el método siguiente denominado método del trapecio como ejemplo de método implícito: \[ \hat{y}_{i+1} =\hat{y}_i+ \frac{h}{2}(f(t_i,\hat{y}_i)+f(t_{i+1},\hat{y}_{i+1})). \]

Proposición.

El método del trapecio es consistente con orden de convergencia \(2\).

Por tanto, aquí se pone de manifiesto lo comentado antes, vemos que el método del trapecio cuya complejidad en su expresión es parecida al método de Euler tiene un orden de convergencia superior a dicho método.

No sólo el método del trapecio es consistente, también es convergente:

Proposición.

Si 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\), el método del trapecio es un método convergente de orden \(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
  • Solve \(\hat{y}_{i+1} =\hat{y}_{i}+\frac{h}{2}(f(t,\hat{y}_{i})+f(t+h,\hat{y}_{i+1}))\) (resolvemos la ecuación anterior con incógnita \(\hat{y}_i\) usando un método numérico para hallar ceros. El valor de \(\hat{y}_i\) será \(\hat{y}(t+ih)\))
  • Set \(t=t+h\). (actualizamos el tiempo \(t\))
• Print \(\mathbf{\hat{y}}\) (Damos el vector de aproximaciones).

Observación.

La ecuación \(\hat{y}_i =\hat{y}_{i-1}+\frac{h}{2}(f(t,\hat{y}_{i-1})+f(t+h,\hat{y}_i))\) a hallar la solución se puede escribir como \(F(\hat{y}_i)=0\), con \(F(x)=\hat{y}_{i-1}+\frac{h}{2}(f(t,\hat{y}_{i-1})+f(t+h,x))-x\).

En caso de querer usar el método de Newton-Raphson, hemos de considerar la sucesión \((\hat{y}_i)_n\) definida por: \[ (\hat{y}_i)_{n+1}=(\hat{y}_i)_{n}-\frac{F((\hat{y}_i)_{n})}{F'((\hat{y}_i)_{n})},\ (\hat{y}_i)_0=\hat{y}_{i-1}. \] El valor de \(F'(x)\) vale: \(F'(x)=\frac{h}{2}\frac{\partial f}{\partial y}(t+h,x)-1.\)

Observación. (continuación)

Así, la sucesión del método de Newton-Raphson será: \[ \begin{align*} (\hat{y}_i)_{n+1} & =(\hat{y}_i)_{n}-\frac{\hat{y}_{i-1}+\frac{h}{2}(f(t,\hat{y}_{i-1})+f(t+h,(\hat{y}_i)_n))-(\hat{y}_i)_n}{\frac{h}{2}\frac{\partial f}{\partial y}(t+h,(\hat{y}_i)_n)-1},\\ & = \frac{\frac{h}{2}\frac{\partial f}{\partial y}(t+h,(\hat{y}_i)_n)(\hat{y}_i)_{n}-\hat{y}_{i-1}-\frac{h}{2}(f(t,\hat{y}_{i-1})+f(t+h,(\hat{y}_i)_n))}{\frac{h}{2}\frac{\partial f}{\partial y}(t+h,(\hat{y}_i)_n)-1},\\ (\hat{y}_i)_0 & =\hat{y}_{i-1}. \end{align*} \]

El estudio de la convergencia del método del trapecio nos permite generalizar el resultado obtenido:

Teorema: relación entre la convergencia y la consistencia de un método de un paso.

Consideremos el siguiente problema de valores iniciales: \[ y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0, \] donde el dominio de definición de \(f\) es de la forma \[ D=\{(t,y)\in\mathbb{R}^2, \ a\leq t\leq b,\ -\infty < y <\infty\}. \] Consideremos el siguiente método de resolución de un paso: \[ \hat{y}_{i+1}=\hat{y}_i +h \varphi(t_i,\hat{y}_i,\hat{y}_{i+1}),\ \hat{y}_0 =y_0. \]

Teorema: relación entre la convergencia y la consistencia de un método de un paso (continuación).

Supongamos que \(\varphi(t,y_i,y_{i+1})\) verifica la condición de Lipschitz en las variables \(y_i\) e \(y_{i+1}\), es decir, existe una constante \(L\) tal que para todo \((t_i,y_i),(t_i,\tilde{y}_i)\in D\), \((t_i,y_{i+1}), (t_i,\tilde{y}_{i+1})\in D\), entonces: \[ |\varphi(t_i,y_i,y_{i+1})-\varphi(t_i,\tilde{y}_i,\tilde{y}_{i+1})|\leq L (|y_i-\tilde{y}_i|+|y_{i+1}-\tilde{y}_{i+1}|). \] Entonces si \(hL<1\), \[ |y(t_i)-\hat{y}_i|\leq \frac{\tau}{2L}\left(\mathrm{e}^{\frac{2(b-a)L}{1-hL}}-1\right), \] donde \(\tau =\max_{0\leq i\leq n}|\tau_i|\).

Además si \(\tau =O(h^k)\), entonces el método converge en orden \(k\).

Observación.

El teorema anterior dice, en pocas palabras, que si la función \(\varphi\) que define el método numérico de un paso es de Lipschitz, entonces si el método es consistente, automáticamente es convergente y además del mismo orden ya que si \(\displaystyle\lim_{h\to 0}\tau =0\), entonces, \[ 0\leq\lim_{h\to 0}|y(t_i)-\hat{y}_i|\leq \frac{\tau}{2L}\left(\mathrm{e}^{\frac{2(b-a)L}{1-hL}}-1\right)=0,\ \Rightarrow \lim_{h\to 0}|y(t_i)-\hat{y}_i|=0. \]

El teorema está escrito de la forma más general posible, suponiendo que el método sea explícito o implícito.

Para calcular el orden de consistencia de un método numérico, necesitamos recordar cuál es la expresión del desarrollo de Taylor de una función de dos variables ya que su cálculo necesita de dicha expresión:

Teorema: desarrollo de Taylor de una función de dos variables.

Sea \(f(t,y)\) una función de dos variables donde sus parciales \(\frac{\partial^{i+j}f}{\partial t^i\partial y^j}(t,y)\) son continuas para todo valor \(i\), \(j\), tal que \(i+j\leq n+1\), con \(n\) un valor entero positivo.

Supongamos que \(f\) está definida en un dominio \(D\) convexo y sea \((t_0,y_0)\in D\).

Entonces para cualquier \((t,y)\in D\), el valor de \(f(t,y)\) puede expresarse de la forma siguiente:

Teorema: desarrollo de Taylor de una función de dos variables (continuación).

\[ f(t,y)=P_n(t,y)+R_n(t,y), \] donde \(P_n(t,y)\) es un polinomio de dos variables de grado menor o igual que \(n\): \[ \begin{align*} P_n(t,y) & = f(t_0,y_0)+\left(\frac{\partial f}{\partial t}(t_0,y_0)(t-t_0)+\frac{\partial f}{\partial y}(t_0,y_0)(y-y_0)\right)\\ & \quad +\frac{1}{2}\left(\frac{\partial^2 f}{\partial t^2}(t_0,y_0)(t-t_0)^2+2\frac{\partial^2 f}{\partial t\partial y}(t_0,y_0)(t-t_0)(y-y_0)\right. \\ &\quad\left. +\frac{\partial^2 f}{\partial y^2}(t_0,y_0)(y-y_0)^2\right)+\cdots \\ & \quad +\frac{1}{n!}\sum_{j=0}^n \binom{n}{j}\frac{\partial^n f}{\partial t^{n-j}\partial y^j}(t_0,y_0)(t-t_0)^{n-j}(y-y_0)^j \end{align*} \]

Teorema: desarrollo de Taylor de una función de dos variables (continuación).

y \(R_n(t,y)\) es el término del error cometido al aproximar \(f\) por el polinomio \(P_n\): \[ R_n(t,y)=\frac{1}{(n+1)!}\sum_{j=0}^{n+1} \binom{n+1}{j}\frac{\partial^{n+1} f}{\partial t^{n-j}\partial y^j}(\xi,\nu)(t-t_0)^{n-j}(y-y_0)^j, \] donde \((\xi,\nu)\) pertenece a un entorno de centro \((t_0,y_0)\) que contiene el valor \((t,y)\).

Observación.

El polinomio \[ P_1(t,y)=f(t_0,y_0)+\left(\frac{\partial f}{\partial t}(t_0,y_0)(t-t_0)+\frac{\partial f}{\partial y}(t_0,y_0)(y-y_0)\right), \] sería la aproximación de primer grado para la función \(f(t,y)\) de error: \[ \begin{align*} R_1(t,y) & = \left(\frac{1}{2}\frac{\partial^2 f}{\partial t^2}(\xi,\nu)(t-t_0)^2+\frac{\partial^2 f}{\partial t\partial y}(\xi,\nu)(t-t_0)(y-y_0)\right. \\ &\quad\left. +\frac{1}{2}\frac{\partial^2 f}{\partial y^2}(\xi,\nu)(y-y_0)^2\right). \end{align*} \]