#!/usr/bin/env python
# coding: utf-8

# # Ecuaciones Diferenciales con Python

# *Esta notebook fue creada originalmente como un blog post por [Raúl E. López Briega](http://relopezbriega.com.ar/) en [Matemáticas, análisis de datos y python](http://relopezbriega.github.io). El contenido esta bajo la licencia BSD.*

# <img alt="Ecuaciones Diferenciales con python" title="Ecuaciones Diferenciales con python" src="http://relopezbriega.github.io/images/DiffEq.png
# " width="500px" height="400px">

# ## Introducción
# 
# Vivimos en un mundo en constante cambio. La posición de la Tierra cambia con el tiempo; la velocidad de un objeto en caída libre cambia con la distancia; el área de un círculo cambia según el tamaño de su *radio*; la trayectoria de un proyectil cambia según la velocidad y el ángulo de disparo.
# Al intentar modelar matemáticamente cualquiera de estos fenómenos, veremos que generalmente adoptan la forma de una o más [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial). 
# 
# **_Nota_**: Antes de continuar leyendo, si no tienen frescos sus conocimientos básicos de [Cálculo integral](https://es.wikipedia.org/wiki/Integraci%C3%B3n) y [Cálculo diferencial](https://es.wikipedia.org/wiki/C%C3%A1lculo_diferencial); les recomiendo que visiten mi anterior artículo de [Introducción al Cálculo](http://relopezbriega.github.io/blog/2015/12/02/introduccion-al-calculo-con-python/).
# 
# ## ¿Qué es una ecuación diferencial?
# 
# Una [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) es una ecuación que involucra una variable dependiente y sus [derivadas](https://es.wikipedia.org/wiki/Derivada) con respecto a una o más variables independientes. Muchas de las leyes de la naturaleza, en [Física](https://es.wikipedia.org/wiki/F%C3%ADsica), [Química](https://es.wikipedia.org/wiki/Qu%C3%ADmica), [Biología](https://es.wikipedia.org/wiki/Biolog%C3%ADa), y [Astronomía](https://es.wikipedia.org/wiki/Astronom%C3%ADa); encuentran la forma más natural de ser expresadas en el lenguaje de las [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial). Estas ecuaciones no sólo tienen aplicaciones en la ciencias físicas, sino que también abundan sus aplicaciones en las ciencias aplicadas como ser [Ingeniería](https://es.wikipedia.org/wiki/Ingenier%C3%ADa), [Finanzas](https://es.wikipedia.org/wiki/Finanzas) y [Economía](https://es.wikipedia.org/wiki/Econom%C3%ADa). 
# 
# Es fácil entender la razón detrás de esta amplia utilidad de las [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial). Si recordamos que $y = f(x)$ es una [función](https://es.wikipedia.org/wiki/Funci%C3%B3n_matem%C3%A1tica), entonces su [derivada](https://es.wikipedia.org/wiki/Derivada) $dy / dx$ puede ser interpretada como el *ritmo de cambio* de $y$ con respecto a $x$. En muchos procesos naturales, las variables involucradas y su *ritmo de cambio* están conectados entre sí por medio de los principios científicos básicos que rigen el proceso. Cuando esta conexión es expresada matemáticamente, el resultado generalmente es una [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial).
# 
# Para ilustrar el caso, veamos un ejemplo. Según la [segunda ley de la dinámica de Newton](https://es.wikipedia.org/wiki/Leyes_de_Newton#Segunda_Ley_de_Newton_o_Ley_Fundamental_de_la_din.C3.A1mica), la [aceleración](https://es.wikipedia.org/wiki/Aceleraci%C3%B3n) $a$ de un cuerpo de [masa](https://es.wikipedia.org/wiki/Masa) $m$ es proporcional a la [fuerza](https://es.wikipedia.org/wiki/Fuerza) total $F$ que actúa sobre él. Si expresamos esto matemáticamente en la forma de una ecuación, entonces tememos que:
# 
# $$F = ma$$
# 
# A primera vista, esta ecuación no parece ser una [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial), pero supongamos que un objeto de [masa](https://es.wikipedia.org/wiki/Masa) $m$ esta en caída libre bajo la influencia de la [fuerza](https://es.wikipedia.org/wiki/Fuerza) de [gravedad](https://es.wikipedia.org/wiki/Gravedad). En este caso, la única [fuerza](https://es.wikipedia.org/wiki/Fuerza) actuando sobre el objeto es $mg$, donde $g$ es la [aceleración](https://es.wikipedia.org/wiki/Aceleraci%C3%B3n) debido a la [gravedad](https://es.wikipedia.org/wiki/Gravedad). Si consideramos a $u$, como la posición del objeto desde una altura determinada; entonces la [velocidad](https://es.wikipedia.org/wiki/Velocidad) de caída del objeto va a ser igual al *ritmo de cambio* de la posición $u$ con respecto al tiempo $t$; es decir que la [velocidad](https://es.wikipedia.org/wiki/Velocidad) es igual a $v = du / dt$, o sea, la [derivada](https://es.wikipedia.org/wiki/Derivada) de la posición del objeto con respecto al tiempo; y como la [aceleración](https://es.wikipedia.org/wiki/Aceleraci%C3%B3n) no es otra cosa que el *ritmo de cambio* de la [velocidad](https://es.wikipedia.org/wiki/Velocidad), entonces podemos definir a la [aceleración](https://es.wikipedia.org/wiki/Aceleraci%C3%B3n) como una segunda [derivada](https://es.wikipedia.org/wiki/Derivada) de la posición del objeto con respecto al tiempo, lo que es igual a decir que $g = d^2u / dt^2$. De esta forma, podemos reescribir la ecuación inicial de la [segunda ley de la dinámica de Newton](https://es.wikipedia.org/wiki/Leyes_de_Newton#Segunda_Ley_de_Newton_o_Ley_Fundamental_de_la_din.C3.A1mica) como la siguiente [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial).
# 
# $$F = m\frac{d^2u}{dt^2}$$
# 
# De esta manera, estamos expresando a la [segunda ley de la dinámica de Newton](https://es.wikipedia.org/wiki/Leyes_de_Newton#Segunda_Ley_de_Newton_o_Ley_Fundamental_de_la_din.C3.A1mica) en función de la posición del objeto. 
# 
# ### Ecuaciones diferenciales ordinarias y parciales
# 
# El caso precedente, es el caso típico de una [Ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria), ya que todas las [derivadas](https://es.wikipedia.org/wiki/Derivada) involucradas son tomadas con respecto a una única y misma variable independiente. Por otro lado, si en la ecuación tenemos [derivadas](https://es.wikipedia.org/wiki/Derivada) de más de una variable independiente, debemos hablar de [Ecuaciones difenciales parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales). Por ejemplo, si $w = f(x, y, z, t)$ es una función de tiempo y 3 coordenadas de un punto en el espacio, entonces las siguientes son sus [Ecuaciones diferenciales parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales):
# 
# $$\frac{\partial^2 w}{\partial x^2} +  \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2} = 0;
# \\
# \\
# a^2\left(\frac{\partial^2 w}{\partial x^2} +  \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2}\right) = \frac{\partial w}{\partial t};
# \\
# \\
# a^2\left(\frac{\partial^2 w}{\partial x^2} +  \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2}\right) = \frac{\partial^2 w}{\partial t^2}.
# $$ 
# 
# En esta caso, el símbolo $\partial$, es la notación matemática para expresar que estamos derivando parcialmente, así por ejemplo, si queremos expresar que vamos a derivar z con respecto a x, escribimos $\partial z / \partial x$.
# 
# Estos ejemplos, tienen una gran aplicación en [Física](https://es.wikipedia.org/wiki/F%C3%ADsica) y se conocen como la [ecuación de Laplace](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_de_Laplace), la [ecuación del calor](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_del_calor) y la [ecuación de onda](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_de_onda) respectivamente. En general, las [Ecuaciones diferenciales parciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_en_derivadas_parciales) surgen en problemas de [campos eléctricos](https://es.wikipedia.org/wiki/Campo_el%C3%A9ctrico), [mecánica de fluidos](https://es.wikipedia.org/wiki/Mec%C3%A1nica_de_fluidos), <a href="https://es.wikipedia.org/wiki/Difusi%C3%B3n_(f%C3%ADsica)">difusión</a> y movimientos de [ondas](https://es.wikipedia.org/wiki/Onda). La teoría de estas ecuaciones es bastante diferente con respecto a las [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) y suelen ser también más complicadas de resolver. En este artículo me voy a enfocar principalmente en las [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria), o [EDOs](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) para abreviar.
# 
# ### Orden de las Ecuaciones diferenciales
# 
# El *orden* de una [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) va a ser igual al *orden* de la mayor [derivada](https://es.wikipedia.org/wiki/Derivada) presente. Así, en nuestro primer ejemplo, la [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) de la [segunda ley de la dinámica de Newton](https://es.wikipedia.org/wiki/Leyes_de_Newton#Segunda_Ley_de_Newton_o_Ley_Fundamental_de_la_din.C3.A1mica) es de segundo orden, ya que nos encontramos ante la segunda derivada de la posición del objeto con respecto al tiempo. 
# La ecuación general de las [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) de grado $n$ es la siguiente:
# 
# $$F\left(x, y, \frac{dy}{dx}, \frac{d^2y}{dx^2}, \dots , \frac{d^ny}{dx^n}\right) = 0 $$
# 
# o utilizando la notación prima para las [derivadas](https://es.wikipedia.org/wiki/Derivada),
# 
# $$F(x, y, y', y'', \dots, y^{(n)}) = 0$$
# 
# La más simple de todas las [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) es la siguiente ecuación de primer orden:
# 
# $$ \frac{dy}{dx} = f(x)$$ 
# 
# y para resolverla simplemente debemos calcular su [integral indefinida](https://es.wikipedia.org/wiki/Integraci%C3%B3n_indefinida):
# 
# $$y = \int f(x) dx + c$$.
# 
# ### Ecuaciones diferenciales separables
# 
# Una *ecuación separable* es una [ecuación diferencial de primer orden](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria_de_primer_orden) en la que la expresión para
# $dx / dy$ se puede factorizar como una función de x multiplicada por una función de y. En otras palabras, puede ser
# escrita en la forma:
# 
# $$\frac{dy}{dx} = f(x)g(y)$$
# 
# El nombre *separable* viene del hecho de que la expresión en el lado derecho se puede "separar" en una función de $x$ y una función de $y$.
# 
# Para resolver este tipo de ecuaciones, podemos reescribirlas en la forma diferencial:
# 
# $$\frac{dy}{g(y)} = f(x)dx$$
# 
# y luego podemos resolver la ecuación original integrando:
# 
# $$\int \frac{dy}{g(y)} = \int f(x) dx + c$$
# 
# Éstas suelen ser las [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) más fáciles de resolver, ya que el problema de resolverlas puede ser reducido a un problema de [integración](https://es.wikipedia.org/wiki/Integraci%C3%B3n); a pesar de que igualmente muchas veces estas integrales pueden ser difíciles de calcular.
# 
# ### Ecuaciones diferenciales lineales
# 
# Uno de los tipos más importantes de [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) son las [Ecuaciones diferenciales lineales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_lineal). Este tipo de ecuaciones son muy comunes en varias ciencias y tienen la ventaja de que pueden llegar a ser resueltas en forma analítica ya que su [ecuación diferencial de primer orden](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria_de_primer_orden) adopta la forma:
# 
# $$\frac{dy}{dx} + P(x)y = Q(x) $$
# 
# donde, $P$ y $Q$ son funciones continuas de $x$ en un determinado intervalo. Para resolver este tipo de ecuaciones de la forma $ y' + P(x)y = Q(x)$, debemos multiplicar los dos lados de la ecuación por el *factor de integración* $e^{\int P(x) dx}$ y luego [integrar](https://es.wikipedia.org/wiki/Integraci%C3%B3n) ambos lados. Así, si por ejemplo quisiéramos resolver la siguiente [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial):
# 
# $$\frac{dy}{dx} + 3x^2y = 6x^2$$
# 
# Primero calculamos el *factor de integración*, $I(x) = e^{\int 3x^2 \ dx}$, lo que es igual a $e^{x^3}$.
# 
# Luego multiplicamos ambos lados de la ecuación por el recién calculado *factor de integración*.
# 
# $$e^{x^3}\frac{dy}{dx} + 3x^2e^{x^3}y = 6x^2e^{x^3}$$
# 
# simplificando, obtenemos:
# 
# $$\frac{d}{dx}(e^{x^3}y) = 6x^2e^{x^3}$$
# 
# Por último, integramos ambos lados de la ecuación:
# 
# $$e^{x^3}y = \int 6x^2e^{x^3} \ dx = 2e^{x^3} + C$$
# 
# y podemos obtener la solución final:
# 
# $$y = Ce^{-x^3} +2 $$
# 
# ### Condición inicial
# 
# Cuando utilizamos [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial), generalmente no estamos interesados en encontrar una familia de soluciones (la solución general), como la que hayamos en el caso precedente, sino que vamos a estar más interesados en una solución que satisface un requerimiento particular.En muchos problemas físicos debemos de encontrar una solución particular que satisface una condición de la forma $y(t_o) = y_0$. Esto se conoce como la **condición inicial**, y el problema de encontrar una solución de la [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) que satisface la *condición inicial* se conoce como el **problema de valor inicial**.
# 
# ## Series de potencias
# 
# Cuando comenzamos a lidiar con las [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial), veremos que existen un gran número de ellas que no pueden ser resueltas en forma analítica utilizando los principios del [Cálculo integral](https://es.wikipedia.org/wiki/Integraci%C3%B3n) y el [Cálculo diferencial](https://es.wikipedia.org/wiki/C%C3%A1lculo_diferencial); pero sin embargo, tal vez podamos encontrar soluciones aproximadas para este tipo de ecuaciones en términos de [Series de potencias](https://es.wikipedia.org/wiki/Serie_de_potencias).
# 
# ### ¿Qué es una serie de potencias?
# 
# Una [Serie de potencias](https://es.wikipedia.org/wiki/Serie_de_potencias) es una [serie](https://es.wikipedia.org/wiki/Serie_matem%C3%A1tica), generalmente infinita, que posee la siguiente forma:
# 
# $$\sum_{n=0}^{\infty} C_nX^n = C_0 + C_1X + C_2X^2 + C_3X^3 + \dots $$
# 
# En dónde $X$ es una variable y las $C_n$ son *constantes* o los *coeficientes* de la serie. Una [Serie de potencias](https://es.wikipedia.org/wiki/Serie_de_potencias) puede [converger](https://es.wikipedia.org/wiki/Serie_convergente) para algunos valores de $X$ y [divergir](https://es.wikipedia.org/wiki/Serie_divergente) para otros valores de $X$. La suma de la [serie](https://es.wikipedia.org/wiki/Serie_matem%C3%A1tica) es una [función](https://es.wikipedia.org/wiki/Funci%C3%B3n_matem%C3%A1tica).
# 
# $$f(x) = C_0 + C_1X + C_2X^2 + C_3X^3 + \dots + C_nX^n + \dots$$
# 
# El *dominio* de esta [función](https://es.wikipedia.org/wiki/Funci%C3%B3n_matem%C3%A1tica) va a estar dado por el conjunto de todos los $X$ para los que la [serie](https://es.wikipedia.org/wiki/Serie_matem%C3%A1tica) [converge](https://es.wikipedia.org/wiki/Serie_convergente).
# 
# ### Series de Taylor
# 
# Las [Series de Taylor](https://es.wikipedia.org/wiki/Serie_de_Taylor) son un caso especial de [Serie de potencias](https://es.wikipedia.org/wiki/Serie_de_potencias) cuyos términos adoptan la forma $(x - a)^n$. Las [Series de Taylor](https://es.wikipedia.org/wiki/Serie_de_Taylor) nos van a permitir aproximar [funciones continuas](https://es.wikipedia.org/wiki/Funci%C3%B3n_continua) que no pueden resolverse en forma analítica y se van a  calcular a partir de las [derivadas](https://es.wikipedia.org/wiki/Derivada) de estas [funciones](https://es.wikipedia.org/wiki/Funci%C3%B3n_matem%C3%A1tica). Su definición matemática es la siguiente:
# 
# $$f(x) = \sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x - a)^n$$
# 
# Lo que es equivalente a decir:
# 
# $$f(x) = f(a) + \frac{f'(a)}{1!}(x -a) + \frac{f''(a)}{2!}(x -a)^2 + \frac{f'''(a)}{3!}(x -a)^3 + \dots$$
# 
# Una de las razones de que las [Series de Taylor](https://es.wikipedia.org/wiki/Serie_de_Taylor) sean importantes es que nos permiten [integrar](https://es.wikipedia.org/wiki/Integraci%C3%B3n) funciones que de otra forma no podíamos manejar. De hecho, a menudo [Newton](https://es.wikipedia.org/wiki/Isaac_Newton) [integraba](https://es.wikipedia.org/wiki/Integraci%C3%B3n) funciones por medio de, en primer lugar, expresarlas como [Series de potencias](https://es.wikipedia.org/wiki/Serie_de_potencias)  y luego [integrando](https://es.wikipedia.org/wiki/Integraci%C3%B3n) la serie término a término. Por ejemplo, la función $f(x) = e^{-x^2}$, no puede ser integrada por los métodos normales del [Cálculo integral](https://es.wikipedia.org/wiki/Integraci%C3%B3n), por lo que debemos recurrir a las [Series de Taylor](https://es.wikipedia.org/wiki/Serie_de_Taylor) para aproximar la solución de su [integral](https://es.wikipedia.org/wiki/Integraci%C3%B3n). Podríamos construir la [Serie de Taylor](https://es.wikipedia.org/wiki/Serie_de_Taylor) de esta función, del siguiente modo:
# 
# $$e^{-x^2} = \sum_{n=0}^{\infty}\frac{(-x^2)^n}{n!} = \sum_{n=0}^{\infty}(-1)^n\frac{x^{2n}}{n!} =
# 1 - \frac{x^2}{1!} + \frac{x^4}{2!} -\frac{x^6}{3!} + \dots$$

# ## Resolviendo Ecuaciones diferenciales con Python
# 
# Mientras que algunos problemas de [Ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) se pueden resolver con métodos analíticos, como hemos mencionado anteriormente, son mucho más comunes los problemas  que no se pueden resolver analíticamente. Por lo tanto, en estos casos debemos recurrir a los métodos numéricos. Es aquí, dónde el poder de las computadoras y en especial, de los paquetes científicos de [Python](https://www.python.org/) como [NumPy](http://www.numpy.org/), [Matplotlib](http://matplotlib.org/), [SymPy](http://www.sympy.org/es/) y [SciPy](http://www.scipy.org/), se vuelven sumamente útiles. Veamos como podemos utilizar la fuerza computacional para resolver [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial).
# 
# ### Soluciones analíticas con Python
# 
# [SymPy](http://www.sympy.org/es/) nos proporciona un solucionador genérico de [Ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria), `sympy.dsolve`, el cual es capaz de encontrar soluciones analíticas a muchas [EDOs](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) elementales. Mientras `sympy.dsolve` se puede utilizar para resolver muchas [EDOs](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) sencillas simbólicamente, como veremos a continuación, debemos tener en cuenta que la mayoría de las [EDOs](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) no se pueden resolver analíticamente. Por ejemplo, retomando el ejemplo que resolvimos analíticamente más arriba, veamos si llegamos al mismo resultado utilizando [SymPy](http://www.sympy.org/es/) para solucionar la siguiente [Ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria):
# 
# $$\frac{dy}{dx} = -3x^2y + 6x^2$$

# In[1]:


# <!-- collapse=True -->
# importando modulos necesarios
get_ipython().run_line_magic('matplotlib', 'inline')

import matplotlib.pyplot as plt
import numpy as np
import sympy 
from scipy import integrate

# imprimir con notación matemática.
sympy.init_printing(use_latex='mathjax') 


# In[2]:


# Resolviendo ecuación diferencial
# defino las incognitas
x = sympy.Symbol('x')
y = sympy.Function('y')

# expreso la ecuacion
f = 6*x**2 - 3*x**2*(y(x))
sympy.Eq(y(x).diff(x), f)


# Aquí primero definimos la incógnitas $x$ utilizando el objeto `Symbol` e $y$, que al ser una función, la definimos con el objeto `Function`, luego expresamos en [Python](https://www.python.org/) la ecuación que define a la función. Ahora solo nos restaría aplicar la función `dsolve` para resolver nuestra [EDO](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria).

# In[3]:


# Resolviendo la ecuación
sympy.dsolve(y(x).diff(x) - f)


# Como podemos ver, arribamos exactamente al mismo resultado. Siguiendo el mismo procedimiento, podemos resolver otras [EDOs](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria), por ejemplo si quisiéramos resolver la siguiente [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial):
# 
# $$\frac{dy}{dx} = \frac{1}{2}(y^2 -1) $$
# 
# que cumpla con la condición inicial $y(0) = 2$, debemos realizar el siguiente procedimiento:

# In[4]:


# definiendo la ecuación
eq = 1.0/2 * (y(x)**2 - 1)

# Condición inicial
ics = {y(0): 2}

# Resolviendo la ecuación
edo_sol = sympy.dsolve(y(x).diff(x) - eq)
edo_sol


# Aquí encontramos la solución general de nuestra [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial), ahora reemplazamos los valores de la condición inicial en nuestra ecuación.

# In[5]:


C_eq = sympy.Eq(edo_sol.lhs.subs(x, 0).subs(ics), edo_sol.rhs.subs(x, 0))
C_eq


# y por último despejamos el valor de la [constante de integración](https://es.wikipedia.org/wiki/Constante_de_integraci%C3%B3n) resolviendo la ecuación.

# In[6]:


sympy.solve(C_eq)


# Como vemos, el valor de la [constante de integración](https://es.wikipedia.org/wiki/Constante_de_integraci%C3%B3n) es 3, por lo tanto nuestra solución para el **problema del valore inicial** es:
# 
# $$y{\left (x \right )} = \frac{3 + e^{x}}{3 - e^{x}}$$
# 
# Para los casos en que no exista una solución analítica, pero sí una solución aproximada por una [Serie de potencias](https://es.wikipedia.org/wiki/Serie_de_potencias), `sympy.dsolve` nos va a devolver la serie como solución. Veamos el caso de la siguiente [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial):
# 
# $$\frac{dy}{dx} = x^2 + y^2 -1$$

# In[7]:


# Solución con serie de potencias
f = y(x)**2 + x**2 -1
sympy.dsolve(y(x).diff(x) - f)


# ### Campos de direcciones
# 
# Los [Campos de direcciones](https://es.wikipedia.org/wiki/Campo_de_direcciones) es una técnica sencilla pero útil para visualizar posibles soluciones a las [ecuaciones diferenciales de primer orden](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria_de_primer_orden). Se compone de líneas cortas que muestran la pendiente de la función incógnita en el plano x-y. Este gráfico se puede producir fácilmente debido a que la pendiente de $y(x)$ en los puntos arbitrarios del plano x-y está dada por la definición misma de la [Ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria):
# 
# $$\frac{dy}{dx} = f(x, y(x))$$
# 
# Es decir, que sólo tenemos que iterar sobre los valores $x$ e $y$ en la grilla de coordenadas de interés y evaluar $f(x, y(x))$ para saber la pendiente de $y(x)$ en ese punto. Cuantos más segmentos de líneas trazamos en un [Campo de dirección](https://es.wikipedia.org/wiki/Campo_de_direcciones), más clara será la imagen. La razón por la cual el gráfico de [Campos de direcciones](https://es.wikipedia.org/wiki/Campo_de_direcciones) es útil, es que las curvas suaves y continuos que son <a href="https://es.wikipedia.org/wiki/Tangente_(geometr%C3%ADa)">tangentes</a> a las líneas de pendiente en cada punto del gráfico, son las posibles soluciones a la [Ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria).
# 
# Por supuesto que calcular las pendientes y dibujar los segmentos de línea para un gran número de puntos a mano sería algo realmente tedioso, pero para eso existen las computadoras y [Python](https://www.python.org/)!. Veamos un ejemplo, supongamos que tenemos la siguiente [Ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria), la cual según lo que vimos más arriba, sabemos que no tiene una solución analítica:
# 
# $$\frac{dy}{dx} = x^2 + y^2 -1$$
# 
# entonces, con la ayuda de [Python](https://www.python.org/), podemos graficar su [Campo de dirección](https://es.wikipedia.org/wiki/Campo_de_direcciones) del siguiente modo:

# In[8]:


# <!-- collapse=True -->
def plot_direction_field(x, y_x, f_xy, x_lim=(-5, 5), y_lim=(-5, 5), ax=None):
    """Esta función dibuja el campo de dirección de una EDO"""
    
    f_np = sympy.lambdify((x, y_x), f_xy, modules='numpy')
    x_vec = np.linspace(x_lim[0], x_lim[1], 20)
    y_vec = np.linspace(y_lim[0], y_lim[1], 20)
    
    if ax is None:
        _, ax = plt.subplots(figsize=(4, 4))
    
    dx = x_vec[1] - x_vec[0]
    dy = y_vec[1] - y_vec[0]
    
    for m, xx in enumerate(x_vec):
        for n, yy in enumerate(y_vec):
            Dy = f_np(xx, yy) * dx
            Dx = 0.8 * dx**2 / np.sqrt(dx**2 + Dy**2)
            Dy = 0.8 * Dy*dy / np.sqrt(dx**2 + Dy**2)
            ax.plot([xx - Dx/2, xx + Dx/2],
                    [yy - Dy/2, yy + Dy/2], 'b', lw=0.5)
    
    ax.axis('tight')
    ax.set_title(r"$%s$" %
                 (sympy.latex(sympy.Eq(y(x).diff(x), f_xy))),
                 fontsize=18)
    
    return ax


# In[9]:


# Defino incognitas
x = sympy.symbols('x')
y = sympy.Function('y')

# Defino la función
f = y(x)**2 + x**2 -1

# grafico de campo de dirección
fig, axes = plt.subplots(1, 1, figsize=(8, 6))
campo_dir = plot_direction_field(x, y(x), f, ax=axes)


# Las líneas de dirección en el gráfico de arriba sugieren cómo las curvas que son soluciones a la [Ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) se van a comportar. Por lo tanto, los [Campos de direcciones](https://es.wikipedia.org/wiki/Campo_de_direcciones) son una herramienta muy útil para visualizar posibles soluciones para [Ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria)  que no se pueden resolver analíticamente. Este gráfico, también nos puede ayudar a determinar el rango de validez de la solución aproximada por la [Serie de potencias](https://es.wikipedia.org/wiki/Serie_de_potencias). Por ejemplo si resolvemos nuestra [EDO](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) para la condición inicial $y(0) = 0$, veamos a que conclusiones arribamos.

# In[10]:


# Condición inicial
ics = {y(0): 0}

# Resolviendo la ecuación diferencial
edo_sol = sympy.dsolve(y(x).diff(x) - f, ics=ics)
edo_sol


# In[11]:


fig, axes = plt.subplots(1, 2, figsize=(10, 5))

# panel izquierdo - solución aproximada por Serie de potencias
plot_direction_field(x, y(x), f, ax=axes[0])
x_vec = np.linspace(-3, 3, 100)
axes[0].plot(x_vec, sympy.lambdify(x, edo_sol.rhs.removeO())(x_vec),
             'b', lw=2)

# panel derecho - Solución por método iterativo
plot_direction_field(x, y(x), f, ax=axes[1])
x_vec = np.linspace(-1, 1, 100)
axes[1].plot(x_vec, sympy.lambdify(x, edo_sol.rhs.removeO())(x_vec),
             'b', lw=2)

# Resolviendo la EDO en forma iterativa 
edo_sol_m = edo_sol_p = edo_sol
dx = 0.125

# x positivos
for x0 in np.arange(1, 2., dx):
    x_vec = np.linspace(x0, x0 + dx, 100)
    ics = {y(x0): edo_sol_p.rhs.removeO().subs(x, x0)}
    edo_sol_p = sympy.dsolve(y(x).diff(x) - f, ics=ics, n=6)
    axes[1].plot(x_vec, sympy.lambdify(x, edo_sol_p.rhs.removeO())(x_vec),
                 'r', lw=2)

# x negativos
for x0 in np.arange(1, 5, dx):
    x_vec = np.linspace(-x0-dx, -x0, 100)
    ics = {y(-x0): edo_sol_m.rhs.removeO().subs(x, -x0)}
    edo_sol_m = sympy.dsolve(y(x).diff(x) - f, ics=ics, n=6)
    axes[1].plot(x_vec, sympy.lambdify(x, edo_sol_m.rhs.removeO())(x_vec),
                 'r', lw=2)


# En el panel de la izquierda podemos ver el gráfico de la solución aproximada por la [Serie de potencias](https://es.wikipedia.org/wiki/Serie_de_potencias); en el panel de la derecha vemos el gráfico al que podemos arribar resolviendo la [EDO](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) para valores incrementales de $x$ para la condición inicial en forma iterativa. Como podemos ver, las solución aproximada se alinea bien con el [campo de direcciones](https://es.wikipedia.org/wiki/Campo_de_direcciones) para los valores de $x$ entre $-1.5$ y $1.5$, luego comienza a desviarse, lo que nos indica que la solución aproximada ya no sería válida. En el panel de la derecha podemos ver una solución que se adapta mucho mejor con el [campo de direcciones](https://es.wikipedia.org/wiki/Campo_de_direcciones).

# ### Transformada de Laplace
# 
# Un método alternativo que podemos utilizar para resolver en forma analítica [Ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) complejas, es utilizar la [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace), que es un tipo particular de [transformada integral](https://es.wikipedia.org/wiki/Transformada_integral). La idea es que podemos utilizar esta técnica para transformar nuestra [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) en algo más simple, resolver esta ecuación más simple y, a continuación, invertir la transformación para recuperar la solución a la [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) original.
# 
# #### ¿Qué es una Transformada de Laplace?
# 
# Para poder comprender la [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace), primero debemos revisar la definición general de la [transformada integral](https://es.wikipedia.org/wiki/Transformada_integral), la cuál adapta la siguiente forma:
# 
# $$T(f(t)) = \int_{\alpha}^{\beta} K (s, t) \ f(t) \ dt = F(s) $$
# 
# En este caso, $f(t)$ es la función que queremos transformar, y $F(s)$ es la función transformada. Los límites de la integración, $\alpha$ y $\beta$, pueden ser cualquier valor entre $-\infty$ y $+\infty$ y $K(s, t)$ es lo que se conoce como el *núcleo* o *kernel* de la transformada, y podemos elegir el *kernel* que nos plazca.  La idea es poder elegir un *kernel* que nos dé la oportunidad de simplificar la [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) con mayor facilidad.
# 
# Si nos restringimos a [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) con *coeficientes* constantes, entonces un *kernel* que resulta realmente útil es $e^{-st}$, ya que al [diferenciar](https://es.wikipedia.org/wiki/C%C3%A1lculo_diferencial) este *kernel* con respecto de $t$, terminamos obteniendo potencias de $s$, que podemos equiparar a los *coeficientes* constantes. De esta forma, podemos arribar a la definición de la [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace):
# 
# $$\mathcal{L}\{f(t)\}=\int_0^{\infty} e^{-st} \ f(t) \ dt$$
# 
# Tengan en cuenta, que además de usar el *kernel* $e^{-st}$, los límites de la integración van desde $0$ a $\infty$, ya que valores negativos de $t$ harían [divergir](https://es.wikipedia.org/wiki/Serie_divergente) la [integral](https://es.wikipedia.org/wiki/Integraci%C3%B3n).
# 
# #### Tabla de transformadas de Laplace
# 
# Calcular la [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) a veces puede resultar en una tarea complicada. Por suerte, podemos recurrir a la siguiente tabla para resolver la [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) para las funciones más comunes:
# 
# <table>
# <thead><tr>
# <th style="text-align:center">Función original</th>
# <th style="text-align:center">Transformada de Laplace</th>
# <th style="text-align:center">Restricciones</th>
# </tr>
# </thead>
# <tbody>
# <tr>
# <td style="text-align:center">$1$</td>
# <td style="text-align:center">$\frac{1}{s}$</td>
# <td style="text-align:center">$s>0$</td>
# </tr>
# <tr>
# <td style="text-align:center">$e^{at}$</td>
# <td style="text-align:center">$\frac{1}{s -a}$</td>
# <td style="text-align:center">$s>a$</td>
# </tr>
# <tr>
# <td style="text-align:center">$t^n$</td>
# <td style="text-align:center">$\frac{n!}{s^{n+1}}$</td>
# <td style="text-align:center">$s>0$, $n$ un entero $> 0$</td>
# </tr>
# <tr>
# <td style="text-align:center">$\cos at$</td>
# <td style="text-align:center">$\frac{s}{s^2 + a^2}$</td>
# <td style="text-align:center">$s>0$</td>
# </tr>
# <tr>
# <td style="text-align:center">$\sin at$</td>
# <td style="text-align:center">$\frac{a}{s^2 + a^2}$</td>
# <td style="text-align:center">$s>0$</td>
# </tr>
# <tr>
# <td style="text-align:center">$\cosh at$</td>
# <td style="text-align:center">$\frac{s}{s^2 - a^2}$</td>
# <td style="text-align:center">$s>|a|$</td>
# </tr>
# <tr>
# <td style="text-align:center">$\sinh at$</td>
# <td style="text-align:center">$\frac{a}{s^2 - a^2}$</td>
# <td style="text-align:center">$s>|a|$</td>
# </tr>
# <tr>
# <td style="text-align:center">$e^{at}\cos bt$</td>
# <td style="text-align:center">$\frac{s- a}{(s^2 - a^2) + b^2}$</td>
# <td style="text-align:center">$s>a$</td>
# </tr>
# <tr>
# <td style="text-align:center">$e^{at}\sin bt$</td>
# <td style="text-align:center">$\frac{b}{(s^2 - a^2) + b^2}$</td>
# <td style="text-align:center">$s>a$</td>
# </tr>
# <tr>
# <td style="text-align:center">$t^n e^{at}$</td>
# <td style="text-align:center">$\frac{n!}{(s-a)^{n+1}}$</td>
# <td style="text-align:center">$s>a$, $n$ un entero $> 0$</td>
# </tr>
# <tr>
# <td style="text-align:center">$f(ct)$</td>
# <td style="text-align:center">$\frac{1}{c}\mathcal{L}\{f(s/c)\}$</td>
# <td style="text-align:center">$c>0$</td>
# </tr>
# </tbody>
# </table>
# 
# 
# #### Propiedades de las transformadas de Laplace
# 
# Las [Transformadas de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) poseen algunas propiedades que también nos van a facilitar el trabajo de resolverlas, algunas de ellas son:
# 
# * ***La [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) es un operador lineal***: Esta propiedad nos dice la [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) de una suma, es igual a la suma de las [Transformadas de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) de cada uno de los términos. Es decir:
# 
# 
# $$\mathcal{L}\{c_1f_1(t) + c_2f_2(t)\}=c_1\mathcal{L}\{f_1(t)\} + c_2\mathcal{L}\{f_2(t)\}$$
# 
# 
# * ***La [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) de la primer derivada***: Esta propiedad nos dice que si $f(t)$ es continua y $f'(t)$ es continua en el intervalo $0 \le t \le \alpha$. Entonces la [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) de la primer derivada es:
# 
# 
# $$\mathcal{L}\{f'(t)\} = s\mathcal{L}\{f(t)\} - f(0)$$
# 
# 
# * ***La [Transformada de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) de derivadas de orden superior***: Esta propiedad es la generalización de la propiedad anterior para derivadas de orden n. Su formula es:
# 
# 
# $$\mathcal{L}\{f^{(n)}(t)\} = s^n\mathcal{L}\{f(t)\} - s^{n-1}f(0) - \dots -f^{(n-1)}(0) \\
# = s^n\mathcal{L}\{f(t)\} - \sum_{i=1}^n s^{n-i}f^{(i-1)}(0)$$
# 
# 
# 
# #### Resolviendo ecuaciones diferenciales con la transformada de Laplace
# 
# La principal ventaja de utilizar [Transformadas de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace) es que cambia la [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) en una [ecuación algebraica](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_algebraica), lo que simplifica el proceso para calcular su solución. La única parte complicada es encontrar las transformaciones y las inversas de las transformaciones de los varios términos de la [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) que queramos resolver. Veamos entonces como [Python](https://www.python.org/) y [SymPy](http://www.sympy.org/es/) nos ayudan a resolver [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) utilizando las [Transformadas de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace). Vamos a intentar resolver la siguiente ecuación:
# 
# $$y'' + 3y' + 2y = 0$$
# 
# con las siguientes condiciones iniciales: $y(0) = 2$ y $y'(0) = -3$

# In[12]:


# Ejemplo de transformada de Laplace
# Defino las incognitas
t = sympy.symbols("t", positive=True)
y = sympy.Function("y")

# Defino la ecuación
edo = y(t).diff(t, t) + 3*y(t).diff(t) + 2*y(t)
sympy.Eq(edo)


# In[13]:


# simbolos adicionales.
s, Y = sympy.symbols("s, Y", real=True)


# In[14]:


# Calculo la transformada de Laplace 
L_edo = sympy.laplace_transform(edo, t, s, noconds=True)
sympy.Eq(L_edo)


# Como podemos ver en este resultado, al aplicar la función `sympy.laplace_transform` sobre la [derivada](https://es.wikipedia.org/wiki/Derivada) de $y(t)$, [SymPy](http://www.sympy.org/es/) nos devuelve un termino con la forma $\mathcal{L}_{t}\left[\frac{d}{d t} y{\left (t \right )}\right]\left(s\right)$. Esto es una complicación si queremos arribar a una [ecuación algebraica](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_algebraica). Por tanto antes de proceder, deberíamos utilizar la siguiente función para resolver este inconveniente.

# In[15]:


def laplace_transform_derivatives(e):
    """
    Evalua las transformadas de Laplace de derivadas de funciones sin evaluar.
    """
    if isinstance(e, sympy.LaplaceTransform):
        if isinstance(e.args[0], sympy.Derivative):
            d, t, s = e.args 
            n = len(d.args) - 1
            return ((s**n) * sympy.LaplaceTransform(d.args[0], t, s) -
                    sum([s**(n-i) * sympy.diff(d.args[0], t, i-1).subs(t, 0)
                         for i in range(1, n+1)]))
        
    if isinstance(e, (sympy.Add, sympy.Mul)):
        t = type(e) 
        return t(*[laplace_transform_derivatives(arg) for arg in e.args])
    
    return e


# In[16]:


# Aplicamos la nueva funcion para evaluar las transformadas de Laplace
# de derivadas
L_edo_2 = laplace_transform_derivatives(L_edo)
sympy.Eq(L_edo_2)


# In[17]:


# reemplazamos la transfomada de Laplace de y(t) por la incognita Y
# para facilitar la lectura de la ecuación.
L_edo_3 = L_edo_2.subs(sympy.laplace_transform(y(t), t, s), Y)
sympy.Eq(L_edo_3)


# In[18]:


# Definimos las condiciones iniciales
ics = {y(0): 2, y(t).diff(t).subs(t, 0): -3}
ics


# In[19]:


# Aplicamos las condiciones iniciales
L_edo_4 = L_edo_3.subs(ics)
L_edo_4


# In[20]:


# Resolvemos la ecuación y arribamos a la Transformada de Laplace
# que es equivalente a nuestra ecuación diferencial
Y_sol = sympy.solve(L_edo_4, Y)
Y_sol


# In[21]:


# Por último, calculamos al inversa de la Transformada de Laplace que 
# obtuvimos arriba, para obtener la solución de nuestra ecuación diferencial.
y_sol = sympy.inverse_laplace_transform(Y_sol[0], s, t)
y_sol


# In[22]:


# Comprobamos la solución.
y_sol.subs(t, 0), sympy.diff(y_sol).subs(t, 0)


# Como podemos ver [Transformadas de Laplace](https://es.wikipedia.org/wiki/Transformada_de_Laplace), pueden ser una buena alternativa para resolver [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) en forma analítica. Pero aún así, siguen existiendo ecuaciones que se resisten a ser resueltas por medios analíticos, para estos casos, debemos recurrir a los métodos numéricos.
# 
# ### Métodos numéricos
# 
# Cuando todo lo demás falla, llegan los métodos numéricos al rescate; siempre vamos a poder calcular una aproximación numérica a una solución de una [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial). De éstos métodos ya no vamos a obtener las formulas elegantes y acabadas que veníamos viendo, sino que vamos a obtener números.
# Existen muchos enfoques para resolver [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) en forma numérica, y generalmente están diseñados para resolver problemas que están formulados como un sistema de [ecuaciones diferenciales de primer orden](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria_de_primer_orden) de la forma estándar:
# 
# $$ \frac{dy}{dx} = f(x, y(x))$$
# 
# donde $y(x)$ es un *[vector](http://es.wikipedia.org/wiki/Vector)* de la funciones incógnitas de $x$.
# 
# La idea básica de muchos de los métodos numéricos para resolver [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) es capturada por el **[Método de Euler](https://es.wikipedia.org/wiki/M%C3%A9todo_de_Euler)**, veamos de que se trata este método.
# 
# #### El Método de Euler
# 
# Para entender este método, revisemos nuevamente la siguiente [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial):
# 
# $$ \frac{dy}{dx} = f(x, y(x))$$
# 
# El [Método de Euler](https://es.wikipedia.org/wiki/M%C3%A9todo_de_Euler) señala que puede que no tengamos la función real que representa la solución a la [Ecuación diferencial](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) precedente. Sin embargo, si poseemos la pendiente de la curva en cualquier punto. Es decir, el *ritmo de cambio* de la curva, que no es otra cosa que su [derivada](https://es.wikipedia.org/wiki/Derivada), la cual podemos utilizar para iterar sobre soluciones en distintos puntos. 
# 
# Supongamos entonces que tenemos el punto $(x_0, y_0)$, que se encuentra en la *curva solución*. Dada la definición de la ecuación que estamos analizando, sabemos que la pendiente de la *curva solución*  en ese punto es $f(x_0, y_0)$. ¿Qué pasaría entonces si quisiéramos encontrar la solución numérica en el punto $(x, y)$ que se encuentra a una corta distancia $h$?. En este caso, podríamos definir a $y$ como $y = y_0 + \Delta y$, es decir, $y$ va ser igual al valor de $y$ en el punto inicial, más el cambio que ocurrió en $y$ al movernos a la distancia $h$. Y como el cambio en $y$ va a estar dado por la pendiente de la curva, que en este caso sabemos que es igual a  $f(x_0, y_0)$, podemos definir la siguiente [relación de recurrencia](https://es.wikipedia.org/wiki/Relaci%C3%B3n_de_recurrencia) que es la base para encontrar las soluciones numéricas con el [Método de Euler](https://es.wikipedia.org/wiki/M%C3%A9todo_de_Euler):
# 
# $$y = y_0 + f(x_0, y_0)h$$
# 
# La cual podemos generalizar para todo $n$ del siguiente forma:
# 
# $$y_n = y_{n -1} + f(x_{n - 1}, y_{n - 1})h$$.
# 
# Este método esta íntimamente relacionado con los [Campos de direcciones](https://es.wikipedia.org/wiki/Campo_de_direcciones) que analizamos más arriba. En general, el [Método de Euler](https://es.wikipedia.org/wiki/M%C3%A9todo_de_Euler) dice que empecemos por el punto dado por el *condición incial* y continuemos en la dirección indicada por el [Campo de direcciones](https://es.wikipedia.org/wiki/Campo_de_direcciones). Luego nos detengamos, miramos a la pendiente en la nueva ubicación, y procedemos en esa dirección.
# 
# #### El Método de Runge-Kutta
# 
# Otro método que podemos utilizar para encontrar soluciones numéricas a [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial), y que incluso es más preciso que el [Método de Euler](https://es.wikipedia.org/wiki/M%C3%A9todo_de_Euler), es el [Método de Runge-Kutta](https://es.wikipedia.org/wiki/M%C3%A9todo_de_Runge-Kutta). En el cual la [relación de recurrencia](https://es.wikipedia.org/wiki/Relaci%C3%B3n_de_recurrencia) va a estar dada por un promedio ponderado de términos de la siguiente manera:
# 
# $$y_{k + 1} = y_k + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$
# 
# en dónde:
# 
# $$k_1 = f(x_k, y_k)h_k, \\
# k_2 = f\left(x_k + \frac{h_k}{2}, y_k + \frac{k_1}{2}\right)h_k, \\
# k_3 = f\left(x_k + \frac{h_k}{2}, y_k + \frac{k_2}{2}\right)h_k, \\
# k_4 = f\left(x_k + h_k, y_k + k_3\right)h_k. \\
# $$
# 
# #### Soluciones numéricas con Python
# 
# Para poder resolver [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial) en forma numérica con [Python](https://www.python.org/), podemos utilizar las herramienta de integración que nos ofrece [SciPy](http://www.scipy.org/), particularmente los dos solucionadores de [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria), `integrate.odeint` y `integrate.ode`. La principal diferencia entre ambos, es que `integrate.ode` es más flexible, ya que nos ofrece la posibilidad de elegir entre distintos *solucionadores*. Veamos algunos ejemplos:
# 
# Comencemos por la siguiente función:
# 
# $$f(x, y(x)) = x + y(x)^2$$

# In[23]:


# Defino la función
f = y(x)**2 + x
f


# In[24]:


# la convierto en una función ejecutable
f_np = sympy.lambdify((y(x), x), f)

# Definimos los valores de la condición inicial y el rango de x sobre los 
# que vamos a iterar para calcular y(x)
y0 = 0
xp = np.linspace(0, 1.9, 100)

# Calculando la solución numerica para los valores de y0 y xp
yp = integrate.odeint(f_np, y0, xp)

# Aplicamos el mismo procedimiento para valores de x negativos
xn = np.linspace(0, -5, 100)
yn = integrate.odeint(f_np, y0, xn)


# Los resultados son dos matrices unidimensionales de [NumPy](http://www.numpy.org/) $yp$ y $yn$, de la misma longitud que las correspondientes matrices de coordenadas $xp$ y $xn$, que contienen las soluciones numéricas de la [ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) para esos puntos específicos. Para visualizar la solución, podemos graficar las matrices $yp$ y $yn$, junto con su [Campo de direcciones](https://es.wikipedia.org/wiki/Campo_de_direcciones).

# In[25]:


fig, axes = plt.subplots(1, 1, figsize=(8, 6))
plot_direction_field(x, y(x), f, ax=axes)
axes.plot(xn, yn, 'b', lw=2)
axes.plot(xp, yp, 'r', lw=2)
plt.show()


# En este ejemplo, solucionamos solo una ecuación. Generalmente, la mayoría de los problemas se presentan en la forma de sistemas de [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria), es decir, que incluyen varias ecuaciones a resolver. Para ver como podemos utilizar a `integrate.odeint` para resolver este tipo de problemas, consideremos el siguiente sistema de [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria), conocido el [atractor de Lorenz](https://es.wikipedia.org/wiki/Atractor_de_Lorenz):
# 
# $$x'(t) = \sigma(y -x), \\
# y'(t) = x(\rho -z)-y, \\
# z'(t) = xy - \beta z
# $$
# 
# Estas ecuaciones son conocidas por sus soluciones caóticas, que dependen sensiblemente de los valores de los parámetros $\sigma$, $\rho$ y $\beta$. Veamos como podemos resolverlas con la ayuda de [Python](https://www.python.org/).

# In[26]:


# Definimos el sistema de ecuaciones
def f(xyz, t, sigma, rho, beta):
    x, y, z = xyz
    return [sigma * (y - x), 
           x * (rho - z) - y,
           x * y - beta * z]

# Asignamos valores a los parámetros
sigma, rho, beta = 8, 28, 8/3.0

# Condición inicial y valores de t sobre los que calcular
xyz0 = [1.0, 1.0, 1.0]
t = np.linspace(0, 25, 10000)

# Resolvemos las ecuaciones
xyz1 = integrate.odeint(f, xyz0, t, args=(sigma, rho, beta))
xyz2 = integrate.odeint(f, xyz0, t, args=(sigma, rho, 0.6*beta))
xyz3 = integrate.odeint(f, xyz0, t, args=(2*sigma, rho, 0.6*beta))


# In[27]:


# Graficamos las soluciones
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig, (ax1,ax2,ax3) = plt.subplots(1, 3, figsize=(12, 4),
                                  subplot_kw={'projection':'3d'})

for ax, xyz, c in [(ax1, xyz1, 'r'), (ax2, xyz2, 'b'), (ax3, xyz3, 'g')]:
    ax.plot(xyz[:,0], xyz[:,1], xyz[:,2], c, alpha=0.5)
    ax.set_xlabel('$x$', fontsize=16)
    ax.set_ylabel('$y$', fontsize=16)
    ax.set_zlabel('$z$', fontsize=16)
    ax.set_xticks([-15, 0, 15])
    ax.set_yticks([-20, 0, 20])
    ax.set_zticks([0, 20, 40])


# Como podemos ver, los solucionadores numéricos que nos ofrece [SciPy](http://www.scipy.org/) son simples de utilizar y pueden simplificar bastante el trabajo de resolver [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial).
# 
# ### Método analítico vs Método numérico
# 
# Al resolver una [ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) en forma analítica, el resultado es una función, $f$, que nos permite calcular la población, $f(t)$, para cualquier valor de $t$. Al resolver una [ecuación diferencial ordinaria](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) en forma numéricamente, se obtienen dos matrices de una dimensión. Podemos pensar a estas matrices como una aproximación discreta de la función continua $f$: *"discreta"*, ya que sólo se define para ciertos valores de $t$, y *"aproximada"*, porque cada valor $F_i$ es sólo una estimación del verdadero valor de $f(t)$. Por tanto, estas son limitaciones de las soluciones numéricas. La principal ventaja es que se puede calcular soluciones numéricas de [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) que no tienen soluciones analíticas, que son la gran mayoría de las [ecuaciones diferenciales ordinarias](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria) no lineales.
# 
# Con esto concluyo este paseo por una de las más fructíferas  ideas de las matemáticas aplicadas, las  [Ecuaciones diferenciales](https://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial); espero que les pueda servir de guía para facilitarles su solución.
# 
# Saludos!
# 
# *Este post fue escrito utilizando Jupyter notebook. Pueden descargar este [notebook](https://github.com/relopezbriega/relopezbriega.github.io/blob/master/downloads/Diff_Eq_python.ipynb) o ver su version estática en [nbviewer](http://nbviewer.ipython.org/github/relopezbriega/relopezbriega.github.io/blob/master/downloads/Diff_Eq_python.ipynb).*