Indice | Previo: Fundamentos.Calculo.CalculoInfinitesimal | Siguiente: Fundamentos.Calculo.CalculoVariacional

4.1.8. Ecuaciones diferenciales

Encontrar la antiderivada de una función (ver la Sección Integrales), se puede formular, de forma general, como el problema de encontrar una función $F(t)$ tal que:

\begin{equation} \label{eq:ecuacion_diferencial_antiderivada} \frac{\mathrm{d}F(t)}{\mathrm{d}t}=f(t) \end{equation}

La solución a este problema es, por definición:

$$ F(t)=\int f(t)\;\mathrm{d}t $$

La integral indefinida en el lado derecho de la anterior ecuación y los métodos numéricos o exactos (analíticos) para obtenerla (no cubiertos en este corto resumen) representan unas de las herramientas matemáticas más útiles de la física.

Pero, existen situaciones en las que el cálculo de una antiderivada no se reduce simplemente a una integral indefinida. Considere por ejemplo el siguiente problema:

\begin{equation} \label{eq:ecuacion_diferencial_simple} \frac{\mathrm{d}^2F(t)}{\mathrm{d}t^2}=-k F(t) \end{equation}

que, en palabras, se formularía como: encontrar la función cuya segunda derivada es proporcional ($k$ se supone constante) al negativo de ella misma.

Ambas, las Ecs. (ecuacion_diferencial_antiderivada) y (ecuacion_diferencial_simple) se conocen como ecuaciones diferenciales.

Las ecuaciones diferenciales se clasifican según:

  • Su orden. El orden de una ecuación diferencial es igual al máximo orden de la derivada de la función objetivo (antiderivada) que aparece en la ecuación. La ecuación diferencial básica (ecuacion_diferencial_antiderivada) es, por ejemplo, de primer orden porque solo involucra la primera derivada de la función $F(t)$. Por su parte, la ecuación diferencial (ecuacion_diferencial_simple) es una ecuación diferencial de segundo orden.

  • Su linealidad. Una ecuación diferencial que solo depende de primeras potencias de la función y sus derivadas se dice que es lineal. En caso contrario tenemos una ecuación diferencial no lineal. Las ecuaciones (ecuacion_diferencial_antiderivada) y (ecuacion_diferencial_simple) son lineales, pero la siguiente ecuación diferencial de primer orden, no lo es:

    $$ \frac{\mathrm{d}F(t)}{\mathrm{d}t}=\frac{h}{F(t)}, $$

    donde $h$ es una constante.

  • El número de variables independientes. Una ecuación diferencial en la que la función depende de una sola variable real se conoce como una ecuación diferencial ordinaria (ODE por la sigla en inglés de ordinary differential equation). Si, por otro lado, la función es un campo escalar o vectorial y la ecuación diferencial se expresa en términos de derivadas parciales (y totales) hablamos de una ecuación diferencial parcial (PDE por sus siglas en inglés).

  • El número de funciones o variables dependientes. Es posible que un problema implique encontrar más de una antiderivada. En ese caso hablamos de un sistema de ecuaciones diferenciales. Un caso común de sistemas de ecuaciones diferenciales se produce cuando queremos encontrar la antiderivada de una función vectorial (cada componente de una función vectorial es una función en sí misma). El caso más importante en la física de un sistema de ecuaciones diferenciales es la ecuación de movimiento de una partícula (que exploraremos a fondo en la Sección Cantidades cinemáticas):

    $$ \frac{\mathrm{d}^2\vec{r}(t)}{\mathrm{d}t^2}=\vec a $$

    Esta ecuación es una forma abrevida de escribir el sistema de ecuaciones diferenciales:

    \begin{eqnarray} \frac{\mathrm{d}^2x(t)}{\mathrm{d}t^2} & = & a_x\nonumber\\ \frac{\mathrm{d}^2y(t)}{\mathrm{d}t^2} & = & a_y\nonumber\\ \frac{\mathrm{d}^2z(t)}{\mathrm{d}t^2} & = & a_z\nonumber\\ \end{eqnarray}

    Como las cantidades $a_x$, $a_y$ y $a_z$ pueden ser a su vez funciones del tiempo, de las funciones $x$, $y$, $z$ y de sus derivadas, se habla, además, de un sistema de ecuaciones diferenciales acopladas.

  • Las condiciones que deben proveerse para resolverla. La solución abstracta de una ecuación diferencial, es decir, el problema de encontrar la antiderivada general, es el equivalente a la integral indefinida. Las integrales definidas, por su lado, equivalen en la teoría de ecuaciones diferenciales a los que se conocen como problemas de valor incial (IVP por el acrónomimo en inglés de initial value problem.) En este tipo de problemas la solución a la ecuación diferencial consiste en encontrar el valor de la función para cualquier valor de la variable independiente una vez se ha provisto el valor de la función (o funciones) y de sus derivadas, en un valor específico o inicial de la variable independiente. Así por ejemplo:

    $$ \frac{\mathrm{d}^2\vec{r}(t)}{\mathrm{d}t^2}=\vec a,\; \;\mathrm{con}\; \vec{r}(0):(0,0,0)\;\mathrm{y}\;\dot{\vec{r}}(0):(1,0,0) $$

    es un IVP.

    Por otro lado un problema de condiciones de frontera (BVP por la sigla en inglés de boundary value problem) es aquel en el que el valor de la función dependiente (no de sus derivadas necesariamente) se provee para varios valores de la variable independiente. Así por ejemplo:

    $$ \frac{\mathrm{d}^2F(t)}{\mathrm{d}t^2}=-F(t),\; \;\mathrm{con}\; F(0)=0\;\mathrm{y}\;F(\pi/2)=1.0, $$

    es un BVP.

Como se intuye fácilmente, la dificultad en la solución a una ecuación diferencial, como sucede también con las ecuaciones algebraicas, puede aumentar con su orden. Sin embargo, usando variables auxiliares, siempre es posible escribir una ecuación diferencial de orden $M$ como un sistema de $M$ ecuaciones diferenciales de primer orden.

Por ejemplo, si en la Ec. (ecuacion_diferencial_simple) llamamos $G(t)\equiv\mathrm{d}F(t)/\mathrm{d}t$, esa ecuación diferencial de segundo orden se puede escribir como el sistema de ecuaciones diferenciales:

\begin{eqnarray} \label{eq:ecuacion_diferencial_simple_reducida} \frac{\mathrm{d}F}{\mathrm{d}t} & = & G\\ \nonumber \frac{\mathrm{d}G}{\mathrm{d}t} & = & -k F\\ \end{eqnarray}

A esta sistema de ecuaciones, lo llamamos el sistema de ecuaciones diferenciales reducido.

La reducción del orden será un método muy utilizado en este libro para abordar la solución a las ecuaciones diferenciales de la mecánica celeste y analítica.

4.1.8.1. Algoritmos para la solución de ODE

La solución aproximada de ecuaciones diferenciales es una de las áreas de mayor interés en el análisis numérico. Sus beneficios prácticos se extienden desde la física teórica y la economía hasta la climatología y la simulación del vuelo de aviones y vehículos espaciales. A lo largo de los últimos 350 años (y en paralelo con la evolución de la mecánica), se han desarrollado métodos numéricos para aproximar la solución de todos los tipos de ecuaciones diferenciales que hemos mencionado hasta aquí.

En este libro, sin embargo, nos concentraremos en la solución de sistemas ecuaciones diferenciales ordinarias con valores iniciales o IVP.

Los métodos numéricos generales, desarrollados para resolver este tipo de problemas (ver Press et al., 2007 para detalles sobre los métodos y algoritmos explícitos), suponen que la ecuación o sistema de ecuaciones diferenciales que queremos resolver puede escribirse como un sistema reducido de ecuaciones diferenciales de primer orden de la forma:

\begin{equation} \label{eq:ecuaciones_reducidas} \{\dot{Y_i} = f_i(\{Y_k\},t)\}_{i=0,1,\ldots,M} \end{equation}

Donde $Y_i$ ($i=0,1,2,\ldots,M-1$) es el conjunto de funciones auxiliares que reemplaza a las funciones dependientes y sus derivadas de orden inferior y $f_i$ son las función que proveen el valor de la primera derivada de la variable auxiliar $Y_i$.

Así por ejemplo, para resolver la ecuación diferencial (ecuacion_diferencial_simple), que ya habíamos reducido como las ecuaciones (ecuacion_diferencial_simple_reducida), las variables auxiliares y sus derivadas serían:

\begin{equation} \label{eq:ecuacion_diferencial_simple_identificacion} \begin{array}{rcl} Y_0=F &,& f_0(t,Y_0,Y_1)=Y_1\\ Y_1=G &,& f_1(t,Y_0,Y_1)=-k Y_0\\ \end{array} \end{equation}

Con esta identificación, el problema original puede escribirse, de forma general como la Ec. (ecuaciones_reducidas).

En este libro usaremos la rutina odeint de la biblioteca científica SciPy (ver Nota abajo) para integrar numéricamente sistemas de ecuaciones diferenciales de la forma reducida. El lector puede leer la documentación completa de odeint para conocer los detalles de su aplicación.

El primer paso para usar odeint es implementar las ecuaciones reducidas como una rutina. En nuestro ejemplo (Ec. ecuacion_diferencial_simple_identificacion) la rutina sería:

In [1]:
def ode_simple(Y,t,k=1):
    f=[0,0]
    f[0]=Y[1]
    f[1]=-k*Y[0]
    return f

Los primeros dos argumentos de esta rutina (ver la Sección Funciones), es decir Y (que contiene una lista de los valores instantáneos de las variables auxiliares $Y_i$) y t (el tiempo en el que las variables auxiliares tienen ese valor) deben estar, estrictamente en ese orden. Otros podrían encontrar más natural poner de primero el tiempo, pero odeint esta diseñado para trabajar con rutinas con este prototipo particular. Además de estos argumentos obligatorios, la rutina puede tener cualquier otro argumento opcional. En este caso aprovechamos esta libertad para proveer el valor de la constante k, que aparece en la ecuación diferencial, y para el cual hemos asumido un valor por defecto k=1 (naturalmente el usuario de la rutina podrá especificar un valor distinto cuando la llame.)

Para resolver este conjunto de ecuaciones diferenciales debemos, además de la rutina anterior, proveer:

  1. Valores específicos para los parámetros de la ecuación diferencial (en este caso la constante $k$),
  2. Una lista de condiciones iniciales, es decir de los valores iniciales de las variables auxiliares $\{Y_i(t=t_0)\}$
  3. un conjunto de valores del tiempo (incluyendo el tiempo inicial $t_0$) para los cuales deseamos predecir el valor de la antiderivada (función o funciones dependientes.)

El siguiente algoritmo prepara estos insumos para odeint en nuestro ejemplo particular:

In [2]:
from numpy import array

k=1.5

Yos=array([1.0,0.0])

ts=array([0.0,1.0,2.0,3.0,4.0,5.0])

Nótese que para las condiciones iniciales y los valores de tiempo (que son aquí arbitrarios, el lector podría escoger unos completamente diferentes) hemos escogido usar arreglos de NumPy (array) en lugar de listas planas (ver la Sección Conjuntos, tuplas y arreglos). Aunque esto no es obligatorio, más adelante hará más fácil la manipulación matemática de estas variables.

Nota: El plural en los algoritmos. Preste atención a la convención que usaremos en lo sucesivo de usar la letra s como sufijo del nombre de algunos arreglos y matrices (p.e. Yos, ts). En lo sucesivo (a no ser que se indique lo contrario) t denotará un valor individual de la variable, pero ts será un arreglo de valores de t.

La solución numérica al cojunto de ecuaciones diferenciales implementados en la rutina ode_simple se obtiene, finalmente, invocando odeint:

In [3]:
from scipy.integrate import odeint
Ys=odeint(ode_simple,Yos,ts,args=(k,))

Solución, Ys =
[[ 1.          0.        ]
 [ 0.33918602 -1.15214115]
 [-0.76990562 -0.78158038]
 [-0.86146852  0.6219388 ]
 [ 0.18550948  1.20348632]
 [ 0.987313    0.1944726 ]]

Las filas de la matriz solución Ys, contienen el valor de las variables auxiliares $\{Y_i\}$ en cada uno de los tiempos provistos. Las columnas, naturalmente, corresponden a los valores instantáneos de cada una de esas variables auxiliares. Así, la componente Ys[0,0] corresponde al valor de $Y_0$ (es decir el valor de la función $F$ de nuestro ejemplo) en $t_0$ (condición inicial).

También es posible extraer tajadas de la matriz. Así, Ys[:,1] (que podría leerse como el segundo valor de cualquier fila o simplemente la columna 1), corresponde al valor de la función auxiliar $G$ en cada uno de los tiempos de integración (recuerde que $G=Y_1$, ver la identificación en la Ec. ecuacion_diferencial_simple_identificacion).

Usando la matriz de solución Ys es posible, finalmente, hacer el gráfico de la función $F(t)$ que estabamos buscando:

In [5]:
%matplotlib inline
In [6]:
import matplotlib.pyplot as plt

#Extraemos los valores de la función F
Fs=Ys[:,0]

plt.figure();
plt.plot(ts,Fs,marker='o',linewidth=0);

#--hide--
plt.xlabel("$t$");
plt.ylabel("$F(t)$");

Figura 4.9. Solución aproximada de la ecuación diferencial $d^2\mathrm{F}/dt^2=-k F$ con k=1.5.

Naturalmente la resolución de este gráfico es bastante pobre porque hemos pedido al algoritmo encontrar únicamente los valores de $F(t)$ en 5 valores del tiempo (arreglo ts.) Si se incrementa el número de componentes de este vector el resultado será mucho más cercano al que esperamos de una función.

Nota: Los algoritmos detrás de odeint. La rutina odeint es un empaque en Python (wrap en inglés) de un complejo y robusto paquete de rutinas conocido como ODEPACK. Desarrollado por el Center for Applied Scientific Computing del Lawrence Livermore National Laboratory, las rutinas de ODEPACK están escritas en lenguaje FORTRAN77 (Python se usa únicamente para pasar los parámetros al paquete y para recuperar las salidas; ese es justamente el sentido del nombre "empaque") y han sido probadas y perfeccionadas durante varias décadas en distintas aplicaciones científicas y de ingeniería (Hindmarsch, 1983).

Existen otras rutinas en el paquete SciPy para resolver ecuaciones diferenciales con condiciones inicales (IVP). Por ejemplo ode y solve_ivp pueden usarse también (esta última es, por ejemplo, la recomendada por los desarrolladores de SciPy). Sin embargo, estas otras rutinas tienen una interface un poco más complicada. Así por ejemplo, para integrar la e.d.m. del ejemplo visto aquí, usando solve_ivp, el código mínimo en Python sería:

from scipy.integrate import solve_ivp
 solucion=solve_ivp(fun=lambda t,Y:ode_simple(Y,t,k),
                    t_span=[ts[0],ts[-1]],y0=Yos,t_eval=ts)

Como puede apreciarse la complejidad del código supera con creces la de aquel que usamos para invocar odeint. A esto se suma el hecho de que la solución, que en el caso de odeint es una matriz Ys fácil de interpretar, en el caso de solve_ivp es en realidad un objeto cuyo atributo solucion.y contiene la solución que buscamos. Y finalmente, pero no menos importante: para el tipo de ecuaciones diferenciales que usaremos en este libro solve_ivp es casi dos veces más lento que odeint. El lector sin embargo puede explorar esas otras alternativas, especialmente si quiere, por ejemplo, comparar distintos métodos de solución (a diferencia de odeint, solve_ivp escoger el método de solución.)

Indice | Previo: Fundamentos.Calculo.CalculoInfinitesimal | Siguiente: Fundamentos.Calculo.CalculoVariacional