Métodos numéricos para resolver ecuaciones diferenciales

Las ecuaciones diferenciales se utilizan para modelar problemas en ciencias e ingeniería que implican el cambio de una variable con respecto a otra. La mayoría de estos problemas requieren la solución de un problema de valor inicial, es decir, la solución de una ecuación diferencial que satisface una condición inicial dada.

En situaciones comunes de la vida real, la ecuación diferencial que modela el problema es demasiado complicada para resolverla exactamente y se adopta uno de dos enfoques para aproximar la solución. El primer enfoque consiste en modificar el problema simplificando la ecuación diferencial a uno que se puede resolver exactamente y luego usar la solución de la ecuación simplificada para aproximar la solución al problema original. El otro enfoque, utiliza métodos para aproximar la solución del problema original. Este es el enfoque que se adopta con mayor frecuencia porque los métodos de aproximación dan resultados más precisos e información de error realista.

Los métodos que se abordarán en esta sección no producen una aproximación continua a la solución del problema de valor inicial. Más bien, las aproximaciones se encuentran en ciertos puntos específicos y, a menudo, igualmente espaciados. Si se necesitan valores intermedios, se utiliza algún método de interpolación.

9.1Método de Euler¶

El método de Euler es la técnica de aproximación más elemental para resolver problemas con valores iniciales. Aunque rara vez se usa en la práctica, la simplicidad de su derivación puede usarse para ilustrar las técnicas involucradas en la construcción de algunas de las técnicas más avanzadas, sin el engorroso álgebra que acompaña a estas construcciones.

El objetivo del método de Euler es obtener aproximaciones al problema de valor inicial:

\frac{dy}{dt} = f(t,y) \,; \qquad a\leq t \leq b \, ; \qquad y(t_0) = y_0

(1)

El método de Euler consiste en ir trazando rectas tangentes sucesivas a partir del valor inicial. Comenzamos linealizando la solución desconocida $y(t)$ de (1) en $t=t_0$ :

L(t) = y_0 + f(t_0, y_0) (t - t_0)

(2)

La gráfica de esta linealización es una línea recta tangente a la gráfica de $y=y(t)$ en el punto $(t_0, y_0)$ y cuya pendiente está dada por $f(t_0, y_0)$ , tal como se puede observar en la Figura Figure 1:

Figure 1:Aproximación mediante el método de Euler

Ahora, considerando un incremento $h$ positivo en el eje $t$ , al reemplazar $t$ por $t_1 = t_0 + h$ en (2) se tiene:

L(t_1) = y_0 + f(t_0, y_0)(t_0 + h - t_0)

(3)

y_1 = y_0 + h f(t_0, y_0)

(4)

Donde $y_1 = L(t_1)$ . En la línea tangente, el punto $(t_1, y_1)$ es una aproximación al punto $(t_1, y(t_1))$ ubicado sobre la curva solución (ver Figura \ref{fig:euler_method}). La exactitud de la aproximación $y_1 \approx y(t_1)$ depende en gran medida del tamaño del incremento $h$ . La recomendación es que debemos elegir este tamaño del paso lo razonablemente pequeño. Ahora repetimos el proceso mediante una segunda línea tangente en $(t_1, y_1)$ . Tomando como punto de base a la aproximación anterior ( $t_1, y_1$ ) podemos obtener una aproximación $y_2 \approx y(t_2)$ correspondiente a dos pasos de longitud $h$ desde $t_0$ , es decir, $t_2 = t_1 + h = t_0 + 2h$ , así:

y(t_2) = y(t_0 + 2h) = y(t_1 + h) \approx y_2 = y_1 + h f(t_1, y_1)

(5)

Si continuamos de esta manera, las sucesivas aproximaciones se pueden expresar de manera general como:

y_{i+1} = y_i + h f(t_i, y_i)

(6)

Donde $t_i = t_0 + ih$ , donde $i=0,1,2,\cdots$ . Observe que cada aproximación se realiza con base en la estimación previa.

Ejemplo Método de Euler

Utilice el método de Euler para obtener una aproximación a $y(0.5)$ , con un paso $h=0.1$ .

\frac{dy}{dt} = 3t + 2 \, ; \qquad y(0)=5

(7)

Solución

Se observa que $f(t,y) = 3t + 2$ , además $t_0=0$ y $y_0 = 5$ . Para la primera aproximación se tiene que:

y_1 = y_0 + h f(t_0, y_0) = 5 + 0.1 (3(0) + 2) = 5.2

(8)

Siendo $t_1 = t_0 + h = 0 + 0.1 = 0.1$ . Para $y_2$ se tiene que:

y_2 = y_1 + h f(t_1, y_1) = 5.2 + 0.1 (3(0.1) + 2) = 5.43 \, ; \qquad t_2 = t_0 + 2h = 0.2

(9)

De manera similar, para las subsiguientes aproximaciones:

y_3 = y_2 + h f(t_2, y_2) = 5.2 + 0.1 (3(0.2) + 2) = 5.69 \, ; \qquad t_3 = t_0 + 3h = 0.3

(10)

y_4 = y_3 + h f(t_3, y_3) = 5.69 + 0.1 (3(0.3) + 2) = 5.98 \, ; \qquad t_4 = t_0 + 4h = 0.4

(11)

y_5 = y_4 + h f(t_4, y_4) = 5.98 + 0.1 (3(0.4) + 2) = 6.3 \, ; \qquad t_5 = t_0 + 5h = 0.5

(12)

Así, la aproximación buscada es $y_{aprox}(0.5) \approx y_5 \approx 6.3$ . Se puede verificar de forma sencilla que la solución exacta para este problema de valor inicial está dada por:

y(t) = \frac{3}{2} t^2 + 2t + 5

(13)

Entonces, el valor real de la función solución en $t=0.5$ es $y(0.5) = 6.375$ . En la Figura Figure 2 se pueden observar la curva solución y la aproximación realizada mediante el método de Euler en el intervalo $0 \leq t \leq 0.5$ . Podemos notar que conforme nos alejamos de la aproximación del valor inicial, la diferencia entre la función solución y la aproximación se hace un poco más evidente.

9.2Métodos de Runge-Kutta¶

Los métodos de Runge-Kutta son una familia de métodos numéricos para aproximar soluciones de ecuaciones diferenciales ordinarias. Son métodos iterativos, lo que significa que se basan en aproximaciones sucesivas para obtener la solución. La idea básica detrás del método de Runge-Kutta es aproximar la solución de una EDO mediante una serie de pasos. En cada paso, se calculan varios valores intermedios que se utilizan para mejorar la aproximación.

Los métodos de Runge-Kutta tienen la forma generalizada:

y_{i+1} = y_{i} + h \, \phi(t_i, y_i, h)

(14)

donde $\phi(t_i, y_i, h)$ se conoce como función incremento, la cual puede interpretarse como una pendiente representativa en el intervalo. La función incremento se escribe en forma general como:

\phi = \left( a_1 k_1 + a_2 k_2 + \cdots + a_n k_n \right)

(15)

donde las $a$ son constantes y las $k$ son:

k_1 = f(t_i, y_i)

(16)

k_2 = f(t_i + p_1 h, y_i + q_{11} k_1 h)

(17)

k_3 = f(t_i + p_2 h, y_i + q_{21} k_1 h + q_{22} k_2 h)

(18)

\vdots \\

(19)

k_n = f(t_i + p_{n-1} h, y_i + q_{n-1,1} k_1 h + q_{n-1,2} k_2 h + \cdots + q_{n-1,n-1} k_{n-1} h)

(20)

donde las $p$ y las $q$ son constantes. Observe que las $k$ son relaciones de recurrencia. Es decir, $k_1$ aparece en la ecuación $k_2$ , la cual aparece en la ecuación $k_3$ , etc. Como cada $k$ es una evaluación funcional, esta recurrencia vuelve eficientes a los métodos RK para cálculos en computadora. Es posible tener varios tipos de métodos de Runge-Kutta empleando diferentes números de términos en la función incremento especificada por $n$ .

9.2.1Método de Runge-Kutta de cuarto orden¶

El más popular de los métodos Runge-Kutta es el de cuarto orden. Hay un número infinito de versiones. La siguiente, es la forma comúnmente usada y, por lo tanto, le llamamos método clásico de Runge-Kutta de cuarto orden (o RK4 en forma abreviada):

y_{i+1} = y_{i} + \frac{h}{6} \left( k_1 + 2k_2 + 2k_3 + k_4 \right)

(21)

donde:

k_1 = f(t_i, y_i)

(22)

k_2 = f\left(t_i + \frac{1}{2}h, y_i + \frac{1}{2}k_1 h \right)

(23)

k_3 = f\left(t_i + \frac{1}{2}h, y_i + \frac{1}{2}k_2 h \right)

(24)

k_4 = f\left(t_i + h, y_i + k_3 h \right)

(25)

Las constantes $k_1$ , $k_2$ , $k_3$ y $k_4$ , representan una pendiente calculada, tal como se muestra en la Figura Figure 3. Cada aproximación es un promedio ponderado de estas pendientes. Observe que $k_2$ se calcula en el medio paso $t_i + 0.5 h$ y utilizando el valor previamente estimado mediante $k_1$ , algo similar ocurre con $k_3$ y $k_4$ .

Figure 3:Pendientes ponderadas del método RK4

Ejemplo. Runge-Kutta de 4to. orden. Utilice un método de Runge-Kutta de cuarto orden para obtener una aproximación al problema de valor inicial en el intervalo $0 \leq t \leq 0.5$ , con un paso $h=0.1$ .

\frac{dy}{dt} = 3t + 2 \, ; \qquad y(0)=5

(26)

Solución

Calculamos las pendientes para la primera aproximación con $t_0 = 0$ y $y_0 = y(t_0) = 5$ :

k_1 = f(t_0, y_0) = 3(0) + 2 = 2

(27)

k_2 = f\left(t_0 + \frac{1}{2}h, y_0 + \frac{1}{2}k_1 h \right) = 3 \rbr{ 0 + 0.5 \rbr{0.1} } + 2 = 0.15 + 2 = 2.15

(28)

k_3 = f\left(t_0 + \frac{1}{2}h, y_0 + \frac{1}{2}k_2 h \right) = 3 \rbr{ 0 + 0.5 \rbr{0.1} } + 2 = 0.15 + 2 = 2.15

(29)

k_4 = f\left(t_0 + h, y_0 + k_3 h \right) = 3 \rbr{ 0 + 0.1 } + 2 = 0.3 + 2 = 2.3

(30)

Entonces, la aproximación $y_1$ está dada por:

y_1 = y_{0} + \frac{h}{6} \left( k_1 + 2k_2 + 2k_3 + k_4 \right) = 5 + \frac{0.1}{6} \rbr{ 2 + 2\rbr{2.15} + 2\rbr{2.15} + 2.3 } = 5.215

(31)

Para la siguiente aproximación $y_2$ , con $t_1 = 0.1$ y $y_1 = 5.215$ , se tiene que:

k_1 = 2.3

(32)

k_2 = 2.45

(33)

k_3 = 2.45

(34)

k_4 = 2.6

(35)

y_2 = 5.46

(36)

Las anteriores y las subsecuentes aproximaciones en el intervalo propuesto se resumen en la siguiente tabla:

$t_i$	$k_1$	$k_2$	$k_3$	$k_4$	$y_i$
0.10	2.0000	2.1500	2.1500	2.3000	5.2150
0.20	2.3000	2.4500	2.4500	2.6000	5.4600
0.30	2.6000	2.7500	2.7500	2.9000	5.7350
0.40	2.9000	3.0500	3.0500	3.2000	6.0400
0.50	3.2000	3.3500	3.3500	3.5000	6.3750

En la Figura Figure 4 se observan las gráficas de las aproximaciones mediante el método de Euler y el de Runge-Kutta de cuarto orden (ambos con $h=0.1$ ). Se puede notar que la aproximación de RK4 prácticamente cae sobre la curva solución real.

9.3Ecuaciones diferenciales ordinarias de segundo orden¶

El problema de valor inicial de segundo orden se puede escribir de forma general como:

\frac{d^2y}{dt^2} = f(t, y, \dot{y})

(37)

y(t_0) = y_0 \, ; \quad \dot{y}(t_0) = u_0

(38)

Este problema se puede escribir como un conjunto de dos ecuaciones diferenciales de primer orden, haciendo la sustitución:

\dot{y} = u

(39)

Entonces:

\dot{u} = f(t,y,u)

(40)

Las ecuaciones (39) y (40) conforman un sistema de dos ecuaciones diferenciales de primer orden. Ahora el problema se centra en determinar las funciones solución $y$ y $u$ sujetas a las condiciones iniciales:

y(t_0) = y_0 \, ; \quad u(t_0) = u_0

(41)

Este sistema de ecuaciones diferenciales se puede resolver numéricamente, aplicando algún método en particular a cada una de las ecuaciones.

9.3.1Método de Euler¶

Así pues, la aproximación mediante el método de Euler para un sistema de EDO de primer orden como el de las ecuaciones (39) y (40), se puede escribir como:

\begin{matrix} y_{i+1} = y_i + h u_i \\[2mm] u_{i+1} = u_i + h f(t_i, y_i, u_i) \end{matrix}

(42)

Ejemplo. Método de Euler, EDO de segundo orden.

Aproxime la solución para:

4\ddot{y} + 4\dot{y} + 17y = 0, \quad y(0)=-1 \, ; \quad \dot{y}(0)=2

(43)

En el intervalo $0 \leq t \leq 1$ , con un paso de $h=0.1$ .

Solución:

La EDO de segundo orden se puede escribir como un sistema de EDO de primer orden mediante la sustitución $\dot{y}=u$ , así:

\dot{y} = u

(44)

\dot{u} = -u - \frac{17}{4} y

(45)

La primera aproximación, con $f(t,y,u)=-u-4.25y$ , $t_0 = 0$ , $y_0 = -1$ y $u_0 = 2$ , está dada por:

y_1 = y_0 + h u_0 = -1 + 0.1(2) = -0.8

(46)

u_1 = u_0 + h f(t_0,y_0,u_0) = 2 + 0.1( -2 - 4.25(-1)) = 2.225

(47)

Para la segunda:

y_2 = y_1 + h u_1 = -0.8 + 0.1(2.225) = -0.5775

(48)

u_2 = u_1 + h f(t_1,y_1,u_1) = 2.225 + 0.1( -2.225 - 4.25(-0.8)) = 2.3425

(49)

Las sucesivas aproximaciones se resumen en la siguiente tabla:

$t_i$	$y_i$	$u_i$
0.00	-1.000000	2.000000
0.10	-0.800000	2.225000
0.20	-0.577500	2.342500
0.30	-0.343250	2.353688
0.40	-0.107881	2.264200
0.50	0.118539	2.083630
0.60	0.326902	1.824888
0.70	0.509390	1.503466
0.80	0.659737	1.136628
0.90	0.773400	0.742577
1.00	0.847658	0.339624

9.4Sistemas de ecuaciones diferenciales de primer orden¶

Considerando que se tiene un sistema de $n$ ecuaciones diferenciales de primer orden de la forma:

\begin{matrix} \displaystyle \frac{dy_1}{dt} = f_1(t, y_1, y_2, \cdots, y_n) \\[2mm] \displaystyle \frac{dy_2}{dt} = f_2(t, y_1, y_2, \cdots, y_n) \\[2mm] \vdots \\[2mm] \displaystyle \frac{dy_n}{dt} = f_n(t, y_1, y_2, \cdots, y_n) \end{matrix}

(50)

Sujetas a las condiciones iniciales:

y_1(t_0)=\alpha_1\,; \quad y_2(t_0) = \alpha_2\,; \quad \cdots \, ; \quad y_n(t_0) = \alpha_n

(51)

El sistema de ecuaciones y las condiciones iniciales se pueden escribir en forma vectorizada como:

\frac{d\vec{y}}{dt} = \vec{f}(t, \vec{y}) \, ; \qquad \vec{y}(t_0) = \vec{y}_0

(52)

Donde:

\vec{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} \, ; \quad \vec{f} = \begin{bmatrix} f_1(t, \vec{y}) \\ f_2(t, \vec{y}) \\ \vdots \\ f_n(t, \vec{y}) \end{bmatrix} \, ; \quad \vec{y}_0 = \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_n \end{bmatrix}

(53)

Observe que la ecuación vectorial (52) es análoga a la ecuación escalar (1). Las expresiones de las secciones previas se pueden utilizar para resolver el sistema de ecuaciones en forma vectorizada, haciendo los ajustes correspondientes en la formulación.

9.4.1Método de Euler¶

Cada aproximación, mediante el método de Euler, del problema de valor de inicial (52) se puede determinar como sigue:

\vec{y}_{i+1} = \vec{y}_i + h \, \vec{f}(t_i, \vec{y}_i)

(54)

donde:

\vec{y}_i = \begin{bmatrix} y_{1_i} \\ y_{2_i} \\ \vdots \\ y_{n_i} \end{bmatrix} \quad ; \quad \vec{f} = \begin{bmatrix} f_1(t, \vec{y}) \\ f_2(t, \vec{y}) \\ \vdots \\ f_n(t, \vec{y}) \end{bmatrix}

(55)

9.4.2Método de Runge-Kutta de cuarto orden¶

Cada aproximación, mediante el método de Runge-Kutta de cuarto orden, del problema de valor de inicial (52) se puede determinar como sigue:

\vec{y}_{i+1} = \vec{y}_i + \frac{h}{6}\left( \vec{k}_1 + 2\vec{k}_2 + 2\vec{k}_3 + \vec{k}_4 \right)

(56)

donde:

\begin{matrix} \vec{k}_1 = \vec{f}(t_i, \vec{y}_i) \\ \vec{k}_2 = \vec{f}\left(t_i + \frac{1}{2}h, \vec{y}_i + \frac{1}{2}\vec{k}_1 h \right) \\ \vec{k}_3 = \vec{f}\left(t_i + \frac{1}{2}h, \vec{y}_i + \frac{1}{2}\vec{k}_2 h \right) \\ \vec{k}_4 = \vec{f}\left(t_i + h, \vec{y}_i + \vec{k}_3 h \right) \end{matrix}

(57)

9.5Sistemas de ecuaciones diferenciales de segundo orden¶

Un sistema de $n$ ecuaciones diferenciales de segundo orden:

\begin{matrix} \displaystyle \ddot{y}_1 = f_1(t, y_1, \dot{y}_1, y_2, \dot{y}_2, \cdots, y_n, \dot{y}_n) \\[2mm] \displaystyle \ddot{y}_2 = f_2(t, y_1, \dot{y}_1, y_2, \dot{y}_2, \cdots, y_n, \dot{y}_n) \\[2mm] \vdots \\[2mm] \displaystyle \ddot{y}_n = f_n(t, y_1, \dot{y}_1, y_2, \dot{y}_2, \cdots, y_n, \dot{y}_n) \\[2mm] \end{matrix}

(58)

Sujeto a las condiciones iniciales:

\begin{matrix} \displaystyle y_1(t_0) = \alpha_1 \\[2mm] \dot{y}_1(t_0) = \alpha_2 \\[2mm] y_2(t_0) = \alpha_3 \\[2mm] \dot{y}_2(t_0) = \alpha_4 \\[2mm] \vdots \\ y_n(t_0) = \alpha_{2n-1} \\[2mm] \dot{y}_n(t_0) = \alpha_{2n} \end{matrix}

(59)

Se puede escribir como un sistema de ecuaciones diferenciales de primer orden, haciendo la sustitución descrita en la sección Section 9.3 para cada una de las ecuaciones, es decir:

\dot{y}_k = u_k

(60)

\dot{u}_k = f_k(t, y_1, u_1, y_2, u_2, \cdots, y_n, u_n)

(61)

Para $k=1, 2, \cdots, n$ . Así pues, el sistema de EDO de segundo orden se puede reescribir como sigue:

\begin{matrix} \displaystyle \dot{y}_1 = u_1 \\ \dot{u}_1 = f_1(t, y_1, u_1, y_2, u_2, \cdots, y_n, u_n) \\ \dot{y}_2 = u_2 \\ \dot{u}_2 = f_2(t, y_1, u_1, y_2, u_2, \cdots, y_n, u_n) \\ \vdots \\ \dot{y}_n = u_n \\ \dot{u}_n = f_n(t, y_1, u_1, y_2, u_2, \cdots, y_n, u_n) \end{matrix}

(62)

Sujeto a las condiciones iniciales:

\begin{matrix} \displaystyle y_1(t_0) = \alpha_1 \\ u_1(t_0) = \alpha_2 \\ y_2(t_0) = \alpha_3 \\ u_2(t_0) = \alpha_4 \\ \vdots \\ y_n(t_0) = \alpha_{2n-1} \\ u_n(t_0) = \alpha_{2n} \end{matrix}

(63)

Un problema de valor de inicial como este se puede resolver de la manera que se describe en la sección Section 9.4.

9.6Problemas¶

Utilice el método de Runge-Kutta de cuarto orden para aproximar la solución del problema de valor inicial de segundo orden. Compare las aproximaciones obtenidas con la solución real.

\ddot{y} - 2\dot{y} + y = te^{t} - t

(64)

y(0) = \dot{y}(0) = 0

(65)

0 \leq t \leq 1 \, ; \quad h=0.1

(66)

\text{Solución real: } y(t) = \frac{1}{6} t^3 e^t - te^t + 2e^t - t - 2

(67)

Utilice el método de RK4 para aproximar la solución al siguiente sistema de ecuaciones diferenciales de primer orden.

\dot{y}_1 = 3y_1 + 2y_2 - \rbr{2t^2 + 1} e^{2t}

(68)

\dot{y}_2 = 4y_1 + y_2 + \rbr{t^2 + 2t -4} e^{2t}

(69)

y_1(0) = 1 ; \quad y_2(0) = 1; \quad 0 \leq t \leq 1 ; \quad h = 0.2

(70)

Compare las aproximaciones obtenidas con las soluciones reales dadas por:

y_1(t) = \frac{1}{3} e^{5t} - \frac{1}{3} e^{-t} + e^{2t}

(71)

y_2(t) = \frac{1}{3} e^{5t} + \frac{2}{3} e^{-t} + t^2 e^{2t}

(72)

Utilice el método de Runge-Kutta de cuarto orden para aproximar la solución al siguiente sistema de ecuaciones diferenciales de primer orden.

\dot{y}_1 = -4y_1 - 2y_2 + \cos t + 4 \sin t

(73)

\dot{y}_2 = 3y_1 + y_2 - 3 \sin t

(74)

y_1(0) = 0 ; \quad y_2(0) = -1; \quad 0 \leq t \leq 2 ; \quad h = 0.1

(75)

Compare las aproximaciones obtenidas con las soluciones reales dadas por:

y_1(t) = 2e^{-t} - 2e^{-2t} + \sin t

(76)

y_2(t) = -3e^{-t} + 2e^{-2t}

(77)

La dinámica de un manipulador simple de un grado de libertad está determinada por la siguiente ecuación diferencial de segundo orden:

m_1 \ddot{q}_1 - m_1 g = \tau_1

(78)

Donde $m_1$ es la masa del manipulador, $g$ la constante de aceleración de la gravedad y $\tau_1$ una fuerza generalizada aplicada al manipulador. Utilice el método de Euler y el método de Runge-Kutta de cuarto orden para aproximar la solución de la ecuación diferencial en el intervalo $0 \leq t \leq 1$ , bajo las siguientes condiciones iniciales y tamaño de paso:

q_1(0) = \dot{q}_1(0) = 0 \, ; \quad h = 0.05

(79)

Considere que $m_1 = 0.5 \text{ kg}$ , $g = 9.81 \text{ m/s}^2$ y $\tau_1 = 5$ N.

9 Métodos numéricos para resolver ecuaciones diferenciales

9.1Método de Euler¶

9.2Métodos de Runge-Kutta¶

9.2.1Método de Runge-Kutta de cuarto orden¶

9.3Ecuaciones diferenciales ordinarias de segundo orden¶

9.3.1Método de Euler¶

9.4Sistemas de ecuaciones diferenciales de primer orden¶

9.4.1Método de Euler¶

9.4.2Método de Runge-Kutta de cuarto orden¶

9.5Sistemas de ecuaciones diferenciales de segundo orden¶

9.6Problemas¶