Ecuaciones diferenciales ordinarias

Introducción

Las ODE son aquellas ecuaciones que responden a la siguiente forma

(1)
\begin{align} \Psi(t, x(t), x'(t), x''(t), \cdots, x^{(n)}(t)) = 0 \end{align}

La anterior es una ecuación diferencial de orden $n$. Cuando en la expresión intervienen sobre la variable independiente funciones no lineales, como el seno o la función exponencial, se dice que es una ecuación diferencial ordinaria no lineal. De lo contrario se la llama lineal.

El objetivo cuando se nos presenta una ecuación diferencial de este tipo es obtener la solución como una expresión de $x(t)$.

Los métodos que se presentan aquí son utilizados para resolver los llamados Initial Value Problem (IVP), donde además de la ecuación diferencial tenemos valores iniciales desde los cuales podemos obtener una aproximación de la solución de la ODE. Muchas veces es conveniente expresar la ecuación 1 de la siguiente manera (ya con los valores iniciales)

(2)
\begin{eqnarray} x^{(n)}(t) &=& \Psi(t, x(t), x'(t), \cdots, x^{(n-1)}(t)) \\ x(t_0) &=& x_0 \\ x'(t_0) &=& x_1 \\ &\vdots& \\ x^{(n-1)}(t_0) &=& x_{n-1} \end{eqnarray}

Ejemplos

  • Ecuación diferencial ordinaria de primer orden lineal
(3)
\begin{align} x'(t) = \lambda x(t) \end{align}
  • Ecuación diferencial ordinaria de segundo orden no lineal
(4)
\begin{align} \ddot{\theta} = - \frac{g}{l} \sin(\theta) \end{align}
  • ODE de orden uno para IVP
(5)
\begin{eqnarray} x'(t) &=& f(t, x(t)) \\ x(0) &=& x_0 \end{eqnarray}

Reducción a un sistema de ecuaciones diferenciales de primer orden

Toda ODE de orden $n$ puede reducirse a un sistema de $n$ ecuaciones diferenciales de primer orden. Para ello se debe asignar una nueva variable a cada una de las derivadas presentes en la expresión.

(6)
\begin{align} \left\{ \begin{array}{lcl} w_1 &=& x(t) \\ w_2 &=& x'(t) \\ &\vdots&\\ w_n &=& x^{(n-1)}(t) \end{array}\right. \end{align}

derivando obtenemos el siguiente sistema

(7)
\begin{align} \left\{ \begin{array}{lcl} w_1' &=& w_2 \\ w_2' &=& w_3 \\ &\vdots&\\ w_n' &=& \Psi(t, x(t), x'(t), \cdots, x^{(n-1)}(t)) \end{array}\right. \end{align}

Los métodos explicados a continuación se basan en ODE de orden 1, para extenderlos es útil reescribir el problema de IVP matricialmente.

(8)
\begin{eqnarray} \mathbf{w}' &=& \mathbf{f}(t, \mathbf{w}) \\ \mathbf{w}(0) &=& \mathbf{w_0} \end{eqnarray}

siendo $\mathbf{w}'$ el vector columna de la derivada de todos los $w_i$, $\mathbf{w_0}$ el vector columna de las condiciones iniciales y $\mathbf{f}(t, \mathbf{w})$ el vector columna de lo que aparece a la derecha del sistema 7.

Método de Euler

Utilizando el Polinomio de Taylor puede obtenerse la siguiente expresión con un orden de error local 2.

(9)
\begin{align} x(t + h) \approx x(t) + h x'(t) \end{align}

Para una ODE de primer orden (ecuación 5) se puede reescribir

(10)
\begin{align} x(t + h) \approx x(t) + h f(t , x(t)) \end{align}

Entonces, eligiendo un paso $h$ podemos obtener una sucesión que representa la solución para el IVP, de la siguiente forma:

(11)
\begin{equation} x_{i+1} = x_i + h f(t_i, x_i) \end{equation}

Ejemplo

Se tiene la siguiente ecuación diferencial de segundo orden, junto con sus valores iniciales:

(12)
\begin{align} \left\{ \begin{array}{l} \theta'' + \theta' + 0.5 \sin(\theta) = 0 \\ \theta(0) = \frac{\pi}{6} \\ \theta'(0) = 0 \end{array}\right. \end{align}

Haciendo el cambio de variables obtenemos un sistema de dos ODE de primer orden (en forma matricial)

(13)
\begin{eqnarray} \mathbf{x} &=& \left( \begin{array}{c} x_1 \\ x_2 \end{array}\right) = \left( \begin{array}{c} \theta \\ \theta' \end{array}\right)\\ \mathbf{x_0} &=& \left( \begin{array}{c} \frac{\pi}{6} \\ 0 \end{array}\right)\\ \mathbf{f}(t,x_1,x_2) &=& \left( \begin{array}{c} x_2 \\ -0.5 \sin(x_1) - x_2 \end{array}\right)\\ \end{eqnarray}

Utilizando el método de Euler (ecuación 11), podemos obtener la siguiente sucesión

(14)
\begin{align} \left( \begin{array}{c} x_1_{i+1} \\ x_2_{i+1} \end{array}\right) = \left( \begin{array}{c} x_1_i \\ x_2_i} \end{array}\right) + h \left( \begin{array}{c} x_2_i \\ -0.5 \sin(x_1_i) - x_2_i \end{array}\right) \end{align}

Para un paso $h = 0.05$ se genera la siguiente tabla de valores:

$i$ $x_{1_i}$ $x_{2_i}$
1 $\frac{\pi}{6}$ 0
2 0.524 -0.013
3 0.523 -0.024
4 0.522 -0.036
10 0.499 -0.099

Entonces, observando los valores que da la columna de $x_{1_i}$ estamos obteniendo una aproximación de la función $\theta$.

Orden del error

El error local de truncamiento (ELT) en el Método de Euler es de orden 2, dado que ese es el orden del error en el Polinomio de Taylor utilizado. Este error es el que se comete al avanzar desde el valor real de la función en un punto al siguiente. El error de truncamiento global, que es el cometido en toda la aproximación de la función, es de orden 1 (siempre es 1 menos que el local).

Referencias

Haga click aquí para ver una linda explicación de cómo se obtiene el Método de Euler.

Método de Heun

Este método es una mejora que se efectúa sobre el Método de Euler. En Euler, se aproxima el siguiente valor de la función conociendo dónde se encuentra ahora y la pendiente en el punto actual, como expresa la ecuación (11). Esto puede llevar a una mala aproximación, pues la pendiente de la función puede cambiar en el intervalo analizado, y Euler arrojaría un valor muy alejado.

Heun propone tomar en consideración tanto la pendiente en el punto actual como la pendiente en el siguiente punto, promediadas, quedando su sucesión de la siguiente manera:

(15)
\begin{align} x_{i+1} = x_i + h \frac{1}{2} ({pend}_{izq} + {pend}_{der}) \end{align}

Pero uno no sabe la pendiente en el siguiente punto dado que uno no conoce exactamente dónde será (si lo supiéramos, ¡no tendría sentido aproximarlo!). Entonces, para aproximar la pendiente "derecha" se utiliza el Método de Euler, como puede verse a continuación.

(16)
\begin{align} \left\{ \begin{array}{lcl} {pend}_{izq} = f(t_i, x_i) \\ {pend}_{der} = f(t_{i+1}, x_i + h f(t_i, x_i)) \end{array}\right. \end{align}

Entonces, se obtiene finalmente la siguiente sucesión:

(17)
\begin{align} x_{i+1} = x_i + h \frac{1}{2} (f(t_i, x_i) + f(t_{i+1}, x_i + h f(t_i, x_i))) \end{align}

Referencias

Haga click aquí para acceder a una muy buena explicación de cómo se obtiene el Método de Heun.

Método de Taylor

Este método utiliza la expansión de Taylor alrededor de un punto y puede alcanzar cualquier orden de error que se desee.

La expansión de Taylor en un punto es:

(18)
\begin{align} x(t + h ) = x(t) + x'(t) h + x''(t) \frac{h^2}{2} + \cdots + x^{(n)}(t) \frac{h^n}{n!} + x^{(n+1)}(\xi) \frac{h^{n+1}}{(n+1)!} \end{align}

Taylor de orden 2

Truncando la expansión de Taylor en el segundo orden se obtiene:

(19)
\begin{align} x(t + h ) = x(t) + x'(t) h + x''(t) \frac{h^2}{2} + x'''(\xi) \frac{h^3}{3!} \end{align}

con $\xi \in [t, t+h]$. De la misma manera que hicimos con Euler para una ODE de orden 1, ahora hay que remplazar $x'(t)$ con $f(t,x(t))$, pero como $f$ es una función de varias variables, cuando aparece $x''(t)$ hay que remplazarlo por la derivada total $f'(t,x(t))$

(20)
\begin{align} x(t + h) = x(t) + h f(t , x(t)) + \frac{h^2}{2} \left( \frac{\partial f (t , x(t))}{\partial t} + \frac{\partial f (t , x(t))}{\partial x} f(t, x(t)) \right) \end{align}

Discretizando podemos obtener la sucesión:

(21)
\begin{align} x_{i+1} = x_i + h f(t_i, x_i) + \frac{h^2}{2} \left( f_t (t_i, x_i) + f_x (t_i, x_i) f(t_i, x_i) \right) \end{align}

Ejemplo (Está mal)

Partiendo del ejemplo de euler tenemos que calcular $\mathbf{f_t}$ y $\mathbf{f_x}$

(22)
\begin{align} \mathbf{f_t} = \left( \begin{array}{c} -0.5 \sin(x_1) - x_2 \\ -0.5 \cos(x_1) x_2 + 0.5 \sin(x_1) + x_2 \end{array}\right) \end{align}

Para calcular $\mathbf{f_x}$, hay que derivar el primer elemento de $f$ por $x_1$ y el segundo por $x_2$ obteniendo:

(23)
\begin{align} \mathbf{f_x} = \left( \begin{array}{c} 0 \\ -1 \end{array}\right) \end{align}

Reemplazando en la ecuación 21 obtenemos:

(24)
\begin{align} \left( \begin{array}{c} x_1_{i+1} \\ x_2_{i+1} \end{array}\right) = \left( \begin{array}{c} x_1_i \\ x_2_i} \end{array}\right) + h \left( \begin{array}{c} x_2_i \\ -0.5 \sin(x_1_i) - x_2_i \end{array}\right) + \frac{h^2}{2} \left( \left( \begin{array}{c} -0.5 \sin(x_1_i) - x_2_i \\ -0.5 \cos(x_1_i) x_2_i + 0.5 \sin(x_1_i) + x_2_i \end{array}\right) + \left( \begin{array}{c} 0 \\ 0.5 \sin(x_1_i) + x_2_i \end{array}\right) \right) \end{align}

Para un paso $h = 0.05$ se genera la siguiente tabla de valores:

$i$ $x_{1_i}$ $x_{2_i}$
1 $\frac{\pi}{6}$ 0
2 0.523 -0.012
3 0.522 -0.023
4 0.519 -0.044
10 0.498 -0.095

Método de Runge-Kutta

Se utilizan las siguientes sucesiones preestablecidas, donde $t_i = i h$:

Orden local 2

(25)
\begin{eqnarray} k_1 &=& h f(t_i, x_i) \\ k_2 &=& h f(t_i + h, x_i + k_1) \\ x_{i+1} &=& x_i + \frac{1}{2} (k_1 + k_2) \end{eqnarray}

Ejemplo

Partiendo del ejemplo de euler tenemos que ver como adaptar el sistema de ecuaciones al algoritmo 25.
La variable $k_1$ se convierte en el vector

(26)
\begin{align} \mathbf{k_1} = \left( \begin{array}{c} h x_2\\ h (-0.5 \sin(x_1) - x_2 ) \end{array}\right) \end{align}

y $k_2$ se obtiene de remplazar en $\mathbf{f}(t, \mathbf{x})$, donde aparezca $t$, poner $t+h$, donde aparezca $x_1$, $x_1 + k_1_1$ siendo $k_1_1$ el primer componente del vector $k_1$ y donde aparezca $x_2$ cambiar por $x_2 + k_1_2$.

(27)
\begin{align} \mathbf{k_2} = \left( \begin{array}{c} h (x_2 + k_1_2)\\ h (-0.5 \sin(x_1 + k_1_1) - (x_2 + k_1_2) ) \end{array}\right) \end{align}

Lo que produce al algoritmo en octave

x1 = pi / 6;
x2 = 0;
h = 0.05;
x = 0;
for i = 1:10
    k11 = h * x2;
    k12 = h * (-0.5 *sin(x1) - x2 );
    k21 = h * (x2 + k12);
    k22 = h * (-0.5 * sin(x1 + k11) - (x2 + k12) );
    x1 = x1 + 0.5 * (k11 + k21);
    x2 = x2 + 0.5 * (k12 + k22);
    printf("Iteracion %d, x1 = %f, x2 = %f\n",i, x1, x2);
end

Para un paso $h = 0.05$ se genera la siguiente tabla de valores:

$i$ $x_{1_i}$ $x_{2_i}$
1 $\frac{\pi}{6}$ 0
2 0.523 -0.012
3 0.522 -0.024
4 0.521 -0.035
10 0.497 -0.097

Orden local 4

(28)
\begin{eqnarray} x_{i+1} &=& x_i + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) \\ k_1 &=& h f(t_i, x_i) \\ k_2 &=& h f(t_i + \frac{h}{2}, x_i + \frac{k_1}{2}) \\ k_3 &=& h f(t_i + \frac{h}{2}, x_i + \frac{k_2}{2}) \\ k_4 &=& h f(t_i + h, x_i + k_3) \end{eqnarray}

Nota: Observar que Runge Kutta de orden 2 es Heun, de la misma manera que Taylor de orden 1 es Euler.

Documentos útiles

Aproximación de soluciones de ecuaciones diferenciales por métodos iterativos

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License