Solución numérica de ecuaciones diferenciales ordinarias

Hecho por: Iván Sánchez & Samuel Ochoa.

En esta libreta se resuelve una serie de ejercicios de ecuaciones diferenciales ordinarias usando los métodos de Euler y de Runge-Kutta de orden 4. Antes de empezar los ejercicios, se presentan las implementaciones en Python de estos métodos:

In [3]:
import numpy as np
import math as m
import matplotlib.pyplot as plt

def euler(f, x0, y0, xf, n):
    """
    Produce la solución aproximada mediante el método 
    de Euler de la EDO y' = f(x, y) con la condición 
    inicial y(x0) = y0, dividiendo el intervalo 
    (x0, xf) en n subintervalos uniformes.
    """
    x    = np.zeros(n)
    y    = np.zeros(n)
    dx   = (xf-x0)/n
    x[0] = x0
    y[0] = y0
    for i in range(n-1):
        x[i+1] = x[i] + dx
        y[i+1] = y[i] + dx * f(x[i], y[i])
    return x, y

def rk4(f, x0, y0, xf, n):
    """
    Produce la solución aproximada mediante el método 
    de Runge-Kutta de orden 4 de la EDO y' = f(x, y) 
    con la condición inicial y(x0) = y0, dividiendo el 
    intervalo (x0, xf) en n subintervalos uniformes.
    """
    x    = np.zeros(n)
    y    = np.zeros(n)
    dx   = (xf-x0)/n
    x[0] = x0
    y[0] = y0
    for i in range(n-1):
        x[i+1] = x[i] + dx
        s1 = f(x[i], y[i])
        s2 = f(x[i] + 0.5*dx, y[i] + 0.5*dx*s1)
        s3 = f(x[i] + 0.5*dx, y[i] + 0.5*dx*s2)
        s4 = f(x[i] + dx, y[i] + dx*s3)
        y[i+1] = y[i] + (dx/6.0)*(s1 + 2*s2 + 2*s3 + s4)
    return x, y

Ejercicio 1

Se sabe que la solución exacta al problema de valor inicial: $$ \left\{ \begin{align*} y' &= \sin(y) \\ y(a) &= y_a \end{align*} \right. $$

está dada por: $$ y(t) = 2\arctan\left(e^{t-a} \tan\left(\frac{y_a}{2}\right)\right) + 2\pi\left\lfloor\frac{y_a + \pi}{2\pi}\right\rfloor $$

donde el símbolo $\lfloor...\rfloor$ denota la función parte entera (piso). Calcule la solución aproximada por Euler del anterior problema de valor inicial:

  • a) tomando $a = 0$, $y_a = 0$
  • b) tomando $a = 0$, $y_a = 100$

En cada caso utilice valores de $\Delta x$ de la forma $0.1\times2^{-k}$ donde $k\in\{0, 1, 2, 3, 4, 5\}$. Realice un gráfico donde compare todas las soluciones numéricas de una misma condición inicial (para cada condición inicial).

1.a.

Tenemos que resolver $y' = \sin(y)$ con la condición inicial $y(0) = 0$. Primero lo intentaremos de forma analítica.

Reescribimos la ecuación y separamos: $$ y' = \sin(y) \\ \frac{dy}{dt} = \sin(y) \\ \frac{dy}{\sin(y)} = dt $$

Integrando ambos lados de la ecuación: $$ \int\!\frac{1}{\sin(y)}\,dy = \int\!dt \\ \int\!\csc(y)\,dy = t \\ \ln\left(csc(y) - \cot(y)\right) = t $$

No resulta aparente una manera de despejar $y$ en términos de $t$, por lo que se demuestra la necesidad de una solución numérica.

Para la solución exacta nos proveen la función $y(t)$: $$ y(t) = 2\arctan\left(e^{t-a} \tan\left(\frac{y_a}{2}\right)\right) + 2\pi\left\lfloor\frac{y_a + \pi}{2\pi}\right\rfloor $$

Implementada en Python sería:

In [4]:
def y1(t, a, ya):
    """Función y(t) generalizada"""
    return (2 * np.arctan(np.exp(t - a) * np.tan(ya / 2.0)) 
           + 2 * np.pi * np.floor((ya + np.pi) / (2 * np.pi)))

def y1a(t):
    """Función y(t) con la condición y(0) = 0"""
    return y1(t, 0, 0)

Tenemos que hallar la solución numérica mediante el método de Euler, por lo que necesitamos definir una función $f(t, y)$ como el lado derecho de $y'$, un intervalo $(x_0, x_f)$ donde se encuentre la solución, y el número divisiones en las que se partirá el intervalo (n).

Vamos a usar el intervalo $(0, 1)$. En lugar del número de divisiones, el enunciado nos pide que utilicemos valores $\Delta x = 0.1 \times 2^{-k}$, donde $k\in\{0,1,2,3,4,5\}$. A partir de $\Delta x$ obtenemos el número de divisiones para cada caso: $$ \Delta x = \frac{x_f - x_0}{n} \Rightarrow n = \frac{x_f - x_0}{\Delta x} $$

A continuación se implementa esta información y se aplica el método de Euler en Python para graficar la solución exacta de la función y sus aproximaciones para cada valor de $\Delta x$:

In [5]:
# Para permitir zoom y desplazamiento.
# Comentar en caso de problemas.
%matplotlib notebook

def graficar1(sol, xmin, xmax, approx, f, x0, y0, xf, 
             divs, fmts, lbls, ejes=False):
    """Función general para graficar PVI con un metodo
    determinado de aproximación numérica.
    """
    x = np.linspace(xmin, xmax)
    y = sol(x)
    fig, ax = plt.subplots()
    if ejes:
        ax.axhline(0.0, color='k')
        ax.axvline(0.0, color='k')
    for i, d in enumerate(divs):
        xx, yy = approx(f, x0, y0, xf, d)
        plt.plot(xx, yy, fmts[i], label=lbls[i])
    ax.plot(x, y, 'k', label='Solución exacta.')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.legend()
    plt.show()

def f1(x, y):
    """Lado derecho de y' = sen(y)"""
    return np.sin(y)

def nd(x0, xf, dx):
    """Número de divisiones dado el intervalo y dx"""
    return int((xf - x0)/dx)

def graficar_1a(xf):
    x0 = 0
    y0 = 0
    dxs = [0.1*2**(-k) for k in range(0, 6)]
    divs = [nd(x0, xf, dx) for dx in dxs]
    fmts = ['r-*','b-v','g-^','y-*','m-v','c-^']
    lbls = [r'$\Delta x$ = ' + str(dx) for dx in dxs]
    graficar1(y1a, x0, xf, euler, f1, x0, y0, xf, 
              divs, fmts, lbls)
    
graficar_1a(1)

Como se puede observar, la solución exacta de este problema es una función constante, y todas las aproximaciones son exactas también. Esto es debido a que no hay variación de un valor a otro que produzca un error.

1.b.

Ahora se nos pide resolver el mismo problema pero con la condición inicial $y(0) = 100$. Como ya tenemos la mayor parte del código implementado, simplemente tenemos que definir la función $y(t)$ particular a la nueva condición inicial y graficar:

In [6]:
def y1b(t):
    return y1(t, 0, 100)

def graficar_1b(xf):
    x0 = 0
    y0 = 100
    dxs = [0.1*2**(-k) for k in range(0, 6)]
    divs = [nd(x0, xf, dx) for dx in dxs]
    fmts = ['r-*','b--','g-.','y:','m--','c-.']
    lbls = [r'$\Delta x$ = ' + str(dx) for dx in dxs]
    graficar1(y1b, x0, xf, euler, f1, x0, y0, xf, 
              divs, fmts, lbls)

graficar_1b(1)

En este caso la solución si varía a lo largo del intervalo, y se puede apreciar que la aproximación por Euler se vuelve más precisa conforme las divisiones del intervalo son más pequeñas o, lo que es lo mismo, conforme el número de divisiones sea mayor. Además, en todos los casos se nota que el método de Euler pierde precisión conforme se aleja del valor inicial.

Ejercicio 2: Euler vs RK4

Considere los siguientes PVI:

  • $y' = 3x(y+4)^2$, $y(0)=0$
  • $y' = e^{x+y}$, $y(0)=0$

Para cada problema de valor inicial:

  • Obtenga la solución analítica usando el método de separación de variables.
  • Utilice los métodos de Euler y RK4 para obtener las soluciones numéricas correspondientes usando $\Delta x = 0.1$
  • Gráfique las soluciones numericas junto con la exacta y compare cuál método funciona mejor en cada caso.

2.a.

Tenemos que resolver el siguiente problema de valor inicial:

$$ \left\{ \begin{align*} y' &= 3x(y+4)^2 \\ y(0) &= 0 \end{align*} \right. $$

Primero resolvemos la ecuación diferencial de forma analítica por separación de variables:

$$ \begin{align*} y' &= 3x(y+4)^2 \\ \frac{dy}{dx} &= 3x(y+4)^2 \\ \frac{dy}{(y+4)^2} &= 3x\,dx \\ \int\!\frac{1}{(y+4)^2}\,dy &= \int\!3x\,dx \end{align*} $$

Para la integral de la izquierda hacemos un cambio de variable:

$$ u = y + 4 \Rightarrow du = dy \\ \\ \begin{align*} \int\!\frac{1}{u^2}\,du &= \int\!3x\,dx \\ \int\!u^{-2}\,du &= \int\!3x\,dx \end{align*} $$

Integrando a ambos lados tenemos la solución general:

$$ \begin{align*} \frac{u^{-1}}{-1} &= 3\cdot\frac{x^2}{2} + C \\ -\frac{1}{u} &= \frac{3}{2}x^2 + C \\ -\frac{1}{y + 4} &= \frac{3}{2}x^2 + C \\ \end{align*} $$

Calculamos la constante de integración dada la condición inicial $y(0) = 0$:

$$ -\frac{1}{0 + 4} = \frac{3}{2}(0)^2 + C \\ \Rightarrow C = -\frac{1}{4} $$

Luego reemplazamos esta constante en la solución general y despejamos la $y$ en términos de $x$ para tener la solución particular de la ecuación diferencial:

$$ \begin{align*} -\frac{1}{y + 4} &= \frac{3}{2}x^2 - \frac{1}{4} \\ \frac{1}{y + 4} &= \frac{1}{4} - \frac{3}{2}x^2 \\ \frac{1}{y + 4} &= \frac{1 - 6x^2}{4} \\ 4 &= (y+4)(1 - 6x^2) \\ \frac{4}{1 - 6x^2} &= y + 4 \\ y &= \frac{4}{1 - 6x^2} - 4 \\ \end{align*} $$

Antes de continuar, como la función es racional, hay que hallar las asíntotas en el eje horizontal para evitarlas. Estas corresponden a los valores para los cuales el denominador es nulo:

$$ 1 - 6x^2 = 0 \Rightarrow x = \pm\sqrt{\frac{1}{6}} \Rightarrow x \approx \pm0,408248290463863 $$

Una vez obtenida esta solución particular exacta del problema, procedemos a obtener la solución numérica por el método de Euler y su error porcentual en cada punto. Nos piden que utilicemos un valor de $\Delta x = 0.1$, por lo que hay que ajustar el número de divisiones (n) de acuerdo a esto. A continuación se presenta la implementación en Python de estos procesos:

In [7]:
import pandas as pd

def f2a(x, y):
    """Lado derecho de y' = f(x, y)"""
    return 3*x * (y + 4)**2

def sol2a(x):
    """Solución particular del problema"""
    return 4/(1 - 6*x**2) - 4

def error_porcentual(x, y, f):
    """Error porcentual del valor aproximado con respecto
    al valor exacto de la función dado por f(x).
    """
    fx = f(x)
    if fx == 0: return 0.0
    return 100.0 * np.absolute((y - f(x)) / f(x))

def generar_tabla(f, x, y, error):
    """Genera una tabla Pandas con los valores aproximados
    (x, y), el valor exacto f(x) y el error relativo a este.
    """
    return pd.DataFrame({
        'x': x,
        'y': y,
        'y(x)': [f(xi) for xi in x],
        'Error relativo (%)': error
    })

def resolver_2a(metodo):
    """Resuelve numéricamente el ejercicio 2.a usando 
    un método determinado.
    """
    x0 = 0.0
    y0 = 0.0
    xf = 0.4  # Para evitar la asíntota en x = 0.408248...
    n  = nd(x0, xf, 0.1)  # Número de divisiones según dx
    x, y = metodo(f2a, x0, y0, xf, n)
    error = [round(error_porcentual(xi, yi, sol2a), 2)
             for xi, yi in zip(x, y)]
    return x, y, error

resultados_2a_euler = resolver_2a(euler)
tabla_2a_euler = generar_tabla(sol2a, *resultados_2a_euler)
tabla_2a_euler
Out[7]:
x y y(x) Error relativo (%)
0 0.0 0.000000 0.000000 0.00
1 0.1 0.000000 0.255319 100.00
2 0.2 0.480000 1.263158 62.00
3 0.3 1.684224 4.695652 64.13

Luego, aplicamos el mismo procedimiento pero con el método de Runge-Kuta de orden 4:

In [8]:
resultados_2a_rk4 = resolver_2a(rk4)
tabla_2a_rk4 = generar_tabla(sol2a, *resultados_2a_rk4)
tabla_2a_rk4
Out[8]:
x y y(x) Error relativo (%)
0 0.0 0.000000 0.000000 0.00
1 0.1 0.255381 0.255319 0.02
2 0.2 1.263307 1.263158 0.01
3 0.3 4.682279 4.695652 0.28

Comparando los valores de las tablas, se nota que el método de Runge-Kuta de orden 4 es mucho más preciso que el método de Euler. A continuación se grafica la solución exacta y ambas aproximaciones numéricas para apreciar mejor la discrepancia:

In [9]:
def graficar2(sol, xmin, xmax, f, x0, y0, xf, n, fmts, ejes=False):
    """Función general para graficar PVI con ambos 
    métodos de aproximación numérica (Euler y RK4).
    """
    x = np.linspace(xmin, xmax)
    y = sol(x)
    fig, ax = plt.subplots()
    if ejes:
        ax.axhline(0.0, color='k')
        ax.axvline(0.0, color='k')
    xe, ye = euler(f, x0, y0, xf, n)
    plt.plot(xe, ye, fmts[0], label='Método de Euler')
    xr, yr = rk4(f, x0, y0, xf, n)
    plt.plot(xr, yr, fmts[1], label='Método de RK4')
    ax.plot(x, y, 'k', label='Solución exacta.')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.legend()
    plt.show()
    
def graficar2a():
    x0 = 0
    y0 = 0
    xf = 0.4
    n = nd(x0, xf, 0.1)
    graficar2(sol2a, x0, xf - 0.05, f2a, x0, y0, xf, n,
             ['r+','b^'])
    
graficar2a()

2.b.

Tenemos que resolver el siguiente problema de valor inicial:

$$ \left\{ \begin{align*} y' &= e^{x+y}\\ y(0) &= 0 \end{align*} \right. $$

Primero resolvemos la ecuación diferencial de forma analítica por separación de variables:

$$ \begin{align*} y' &= e^{x+y} \\ \frac{dy}{dx} &= e^{x+y} \\ \frac{dy}{dx} &= e^{x}\,e^{y} \\ \frac{dy}{e^{y}} &= e^{x}\,dx \\ \int\!{e^{-y}}\,dy &= \int\!e^{x}\,dx \end{align*} $$

Integrando a ambos lados tenemos la solución general:

$$ \begin{align*} \frac{e^{-y}}{-1} &= {e^{x}} + C \\ -e^{-y} &= {e^{x}} + C \\ \end{align*} $$

Calculamos la constante de integración dada la condición inicial $y(0) = 0$:

$$ -e^{-0} = e^{0} + C \\ -1 = 1 + C \\ \Rightarrow C = -2 $$

Luego reemplazamos esta constante en la solución general y despejamos la $y$ en términos de $x$ para tener la solución particular de la ecuación diferencial:

$$ \begin{align*} -e^{-y} &= {e^{x}} - 2 \\ e^{-y} &= -{e^{x}} + 2 \\ ln(e^{-y}) &= ln(2 -{e^{x}})\\ {-y} &= ln(2 -{e^{x}}) \\ {y} &= -ln(2 -{e^{x}}) \\ \end{align*} $$

Antes de continuar, hay que tomar en cuenta que la función logaritmo neperiano solo está definida para valores en los que su argumento es mayor a cero:

$$ 2 -{e^{x}} > 0 \Rightarrow {e^{x}} < 2 \Rightarrow x < ln(2) \Rightarrow x < 0,6931471805599453 \text{ (aprox)} $$

Entonces, la función solución está definida para valores menores a $ln(2)$, por lo que el intervalo no puede pasar dicho valor. Con esta información, procedemos a obtener la solución numérica por el método de Euler y su error porcentual en cada punto. Nos piden que utilicemos un valor de $\Delta x = 0.1$, por lo que hay que ajustar el número de divisiones (n) de acuerdo a esto. A continuación se presenta la implementación en Python de estos procesos:

In [10]:
def f2b(x, y):
    """Lado derecho de y' = f(x, y)"""
    return np.exp(x + y)

def sol2b(x):
    """Solución particular del problema"""
    return -np.log(2 - np.e**x)

def resolver_2b(metodo):
    """Resuelve numéricamente el ejercicio 2.a usando 
    un método determinado.
    """
    x0 = 0.0
    y0 = 0.0
    xf = 0.6  # x tiene que ser menor a 0.693147...
    n  = nd(x0, xf, 0.1)  # Número de divisiones según dx
    x, y = metodo(f2b, x0, y0, xf, n)
    error = [round(error_porcentual(xi, yi, sol2b), 2)
             for xi, yi in zip(x, y)]
    return x, y, error

resultados_2b_euler = resolver_2b(euler)
tabla_2b_euler = generar_tabla(sol2b, *resultados_2b_euler)
tabla_2b_euler
Out[10]:
x y y(x) Error relativo (%)
0 0.00 0.000000 -0.000000 0.00
1 0.12 0.120000 0.136389 12.02
2 0.24 0.272550 0.316423 13.87
3 0.36 0.472895 0.567977 16.74
4 0.48 0.748891 0.957307 21.77

Luego, aplicamos el mismo procedimiento pero con el método de Runge-Kuta de orden 4:

In [11]:
resultados_2b_rk4 = resolver_2b(rk4)
tabla_2b_rk4 = generar_tabla(sol2b, *resultados_2b_rk4)
tabla_2b_rk4
Out[11]:
x y y(x) Error relativo (%)
0 0.00 0.000000 -0.000000 0.00
1 0.12 0.136390 0.136389 0.00
2 0.24 0.316427 0.316423 0.00
3 0.36 0.567989 0.567977 0.00
4 0.48 0.957357 0.957307 0.01

Comparando los valores de las tablas, se nota que el método de Runge-Kuta de orden 4 es mucho más preciso que el método de Euler. Además, el método de Euler se vuelve cada vez menos preciso conforme se aleja del punto inicial. A continuación se grafica la solución exacta y ambas aproximaciones numéricas para apreciar mejor la discrepancia:

In [12]:
def graficar2b():
    x0 = 0
    y0 = 0
    xf = 0.6
    n = nd(x0, xf, 0.1)
    graficar2(sol2b, x0, xf, f2b, x0, y0, xf, n,
             ['r+','b^'])
    
graficar2b()
In [ ]: