Practica_EDOs_2011_12

Introducción a las Ecuaciones Diferenciales Ordinarias

1. Motivación y nociones básicas

En las Ciencias y la Ingeniería se desarrollan modelos matemáticos para comprender mejor los fenómenos físicos. Muchos de estos modelos provienen de fenómenos de la vida real en los que una magnitud cambia de estado con el tiempo y proporcionan una ecuación que contiene alguna(s) derivada(s) de una función desconocida o incógnita. Esta ecuación es una ecuación diferencial. Si la ecuación contiene derivadas parciales de la función incógnita, la ecuación se denomina ecuación en derivadas parciales (EDP). Si en la ecuación solo aparecen derivadas de una función desconocida de una variable, la ecuación se denomina ecuación diferencial ordinaria. Nos vamos a centrar en este último tipo de ecuaciones.

Más precisamente, una  ecuación diferencial ordinaria (EDO, en adelante) es una ecuación que relaciona una función desconocida o incógnita, y(x), con su variable independiente, x, y con una o más de sus derivadas hasta un cierto orden. En muchos casos prácticos, la variable independiente (v.i.) es el tiempo y se denota por t. Nosotros utilizaremos ambas notaciones para la v.i.

Se llama orden de una EDO al orden de la mayor derivada de la función incógnita que aparece en la ecuación. 

 Ejemplo 1.  Las  ecuaciones diferenciales  

son, respectivamente, de tercer, primer y segundo orden.

En adelante, denotaremos y=y(x), y'=y'(x), y''=y''(x),...

Resolver una EDO consiste en encontrar, si es posible, una función (o funciones) y=y(x) que verifique(n) la EDO para cualquier x en un cierto intervalo real I. Se dirá que dicha función (o funciones) es (son) solución (soluciones) de la EDO en el intervalo I. La resolución de una EDO conlleva el cálculo de alguna integral. Por ello, a la resolución de una EDO también se le llama integración de una EDO. Resaltar que hay EDOs que no tienen solución, como la EDO

2. Un ejemplo sencillo

Ejemplo 2. Resolver (o integrar) la EDO de primer orden 

Se trata de hallar una función y=y(x) tal que y'(x)  = 2x (o bien \displaystyle \frac{dy}{dx} = 2x) para cualquier x en algún intervalo real I. Es claro que la función y(x) = x^2 es una solución de la ecuación (ed) para cualquier x real. Pero también es solución de (ed) cualquier  función de la forma

De hecho, para resolver (ed), basta con integrar

Esta fórmula, que contiene a las infinitas soluciones de (ed), se llama solución general de (ed). Para cada valor de la constante C se obtiene una solución diferente. Cada una de ellas se llama solución particular (o curva integral) de (ed).

Este ejemplo es un caso particular del tipo de EDOs que pasamos a estudiar.

3. La EDO y'=f(x)

Consideramos la EDO de primer orden

donde f es una función continua en su dominio de definición D = Dom \, (f). La fórmula que proporciona el cálculo

integral de las infinitas primitivas de f, 

Las primitivas conocidas (que se expresan con un número finito de funciones elementales) se obtienen con cualquier programa razonable de cálculo simbólico. En particular, también con Sage, como veremos.

4. La EDO de primer orden y' = f(x,y)

Las EDOs de primer orden son las más sencillas de estudiar.  Las podemos expresar en su forma implícita:

o en su forma normal (o explícita):

Ejemplo 3. La EDO 1 + y^2 + x y y' = 0 se puede escribir en forma normal (sin más que despejar y' de la expresión dada):   

Nos centraremos en las EDOs de primer orden escritas en forma normal.

Una solución particular (o curva integral) de la EDO  

Ejercicio: Sin hacer uso del ordenador, compruébese que \varphi_1(x)=e^{-2x}+1-2x+2x^2 es una solución en toda la recta real de la EDO

4.1 Obtención con Sage de la solución general de una EDO de primer orden

Los métodos de resolución de EDOs son muy limitados y, en general, no cabe esperar resolver de forma explícita una ecuación diferencial ordinaria cualquiera dada. Se conoce solo una colección finita de tipos de EDOs con método explícito de resolución, como las ecuaciones diferenciales de primer orden lineales (homogéneas y no homogéneas), las de Bernoulli, las exactas, las reducibles a exactas con factor integrante, las homogéneas y algunas reducibles a ellas... La mayoría de ellos están implementados en los programas usuales de cálculo simbólico.

Para resolver con Sage una EDO de primer orden, tenemos que:

• declarar la variable independiente x;
• declarar la función incógnita como una función que depende de x;
• almacenar la ecuación diferencial en una variable (por ejemplo, de), utilizando la orden diff(y,x) para la derivada de y respecto de x;
• resolver la ecuación diferencial utilizando la orden desolve.  La sintaxis de esta orden es:

Observaciones:

a) Al almacenar la ecuación diferencial, póngase cuidado en no confundir la asignación "=" con el signo de igualdad "==".

b) Se pueden añadir otros atributos a la orden desolve.

c) Hay otros métodos de resolución de EDOs con Sage (como desolve_Laplace,...),  que no estudiaremos aquí.

4.2. Algunos ejemplos de resolución de EDOs de primer orden

Ejemplo 4. Comenzamos resolviendo la EDO 

También podemos resolver la ecuación anterior almacenándola en la forma

1. Asimismo, podría usarse directamente la orden desolve sin almacenar previamente la EDO en una variable. Aquí no añadimos el atributo show_method, que era opcional.

Ejemplo 5. Es claro que cualquier función de la forma \frac{1}{2} e^{2x} +C es solución de la EDO y' = e^{2x} (de hecho, es su solución general). Lo comprobamos con Sage: 

Con desolve podemos resolver ecuaciones diferenciales ordinarias de primer orden, que no sean necesariamente de la forma y'=f(x).

Ejemplo 6. Resolvamos las siguientes EDOs de primer orden:

1) y'= y \qquad 2) y'= 7y\qquad 

(ejemplos de ecuaciones de primer orden lineales homogéneas).

3) 2y'+y=3 \qquad 4) y'-2xy=x

(ejemplos de ecuaciones de primer orden lineales no homogéneas)

5) y' + y = x y^2 (ejemplo de ecuación de Bernoulli)

6) yy' + 3x=0 (ejemplo de ecuación de variables separables)

Si la EDO está en forma implícita, como es este caso, al almacenar la ecuación diferencial no es necesario explicitar en el segundo miembro que está igualada a 0. En efecto, podíamos haber hecho:

7) y^2y' + y =3 (ejemplo de ecuación autónoma: la v.i. no aparece de manera explícita).

8) x+y+ (x-y) y'=0 (ejemplo de ecuación exacta: de la forma M(x,y) + N(x,y) y' = 0.)

¡Ojo!

• Obsérvese que en los tres últimos ejemplos, la familia solución viene expresada en forma implícita.
• EDOs difíciles pueden producir error, como vemos al intentar resolver con Sage la EDO  
  \sqrt{y}y' + e^y +\cos x - {\rm sen} (x+y)=0.

• De hecho, en algunos casos, Sage no es capaz de distinguir si una EDO tiene solución o no, como es el caso de la EDO
  1 + (y')^2 = 0,
  que no tiene solución. Veámoslo con Sage.

Ejemplo 7. Resolver el problema de valor inicial

(PVI) \quad   \left\{ \begin{array}{l} 2y' + y =3,  \\ y(1)=3. \end{array} \right.  

De hecho, puede demostrarse que cualquier (PVI) de la forma

(PVI) \quad   \left\{ \begin{array}{l} y'  =y^{\alpha} \\ y(x_0)=0, \end{array} \right.  

con \alpha \in (0,1), tiene infinitas soluciones (una de ellas la solución nula). No lo vemos aquí.

4.4. Dibujando curvas integrales

(nótese que hemos tenido que declarar la variable c). A continuación, como vimos en la Práctica 1, con el método .plot podemos representar varias curvas integrales, por ejemplo,

en el intervalo [-4,4], para los valores c=1,2,3,4,5 (por ejemplo) de la constante c.

• Podemos dibujar soluciones particulares mediante el método append. Veámoslo con el siguiente ejemplo.

Ejemplo 9. Dibujamos algunas soluciones particulares de la EDO

Para apreciar mejor el comportamiento de las curvas integrales (o soluciones particulares) cerca del origen, podemos hacer un zoom y representar, por ejemplo, soluciones particulares con la variable dependiente y entre -10 y 10:

O con un zoom mayor, por ejemplo, con y \in [-2,2]:

4.5. Campo de direcciones

En general, las soluciones de una EDO de primer orden no se pueden hallar de forma explícita, ni siquiera en el caso lineal (pues las integrales involucradas pueden tener primitivas no elementales). Por ello, se suelen utilizar esquemas gráficos, para visualizar su comportamiento, y métodos numéricos, para aproximar sus valores.

Interpretación geométrica de la solución de un (PVI)

Supongamos dado un 

  (PVI) \quad   \left\{ \begin{array}{l} y'=f(x,y)  \\ y(x_0)=y_0, \end{array} \right.  

con solución única \varphi(x) en un intervalo real I que contine a x_0 en su interior. En particular,  y= \varphi(x) verifica

es decir, la curva y= \varphi(x) pasa por el punto (x_0,y_0) y la pendiente de la recta tangente a y= \varphi(x) que pasa por dicho punto (x_0,y_0) tiene pendiente  f(x_0,y_0). De este modo, aunque no conozcamos la expresión explícita de la solución \varphi(x) podemos hacer un esbozo de ella "cerca del punto (x_0,y_0)",  dibujando su aproximación lineal mediante un segmento que pase por (x_0,y_0) con pendiente \varphi'(x_0)=f(x_0,y_0).

Si esto lo hacemos con cada punto (x_0,y_0) de un recinto dado del plano, obtenemos el campo direccional o campo de direcciones.

Representación gráfica del campo de direcciones

Los programas usuales de cálculo numérico o simbólico suelen implementar el dibujo de campos direccionales. Sage dibuja el campo de direcciones con la orden plot_vector_field. Veamos su sintaxis con la EDO del ejemplo anterior:

Nota importante: Si la pendiente de la recta tangente a la curva y= \varphi(x) que pasa por el punto (x,y) es f(x, y), un vector director de dicha recta es (1, f(x,y)) o también  (\lambda, \lambda f(x,y)) para cualquier \lambda real no nulo. Como en el ejemplo que vamos a tratar es f(x,y) = 2y/x + x^2, se toma en el punto (x,y) el vector director (x, xf(x,y)).

Podemos también representar gráficamente las curvas integrales junto con el campo de direcciones. Lo hacemos para el mismo ejemplo:

Representemos ahora en el rectángulo [0,1]\times [0,1] curvas integrales y el campo de direcciones de la EDO

Finalmente, para la misma EDO, 2y' + y -3 = 0, vamos a representar en el rectángulo [0,8] \times [-2,6] el campo direccional y las curvas integrales (o soluciones particulares) que pasan por los puntos (0,-2), (0,0), (0,2), (0,4) y (0,6).

4.6. Métodos numéricos de resolución de EDOs de primer orden: Métodos de Euler y de Heun

Nuevamente, supongamos dado un 

  (PVI) \quad   \left\{ \begin{array}{l} y'=f(x,y)  \\ y(x_0)=y_0, \end{array} \right.  

con solución única \varphi(x) en un intervalo real I que contine a x_0 en su interior. Resolver numéricamente el (PVI) dado consiste en considerar un conjunto finito 

de puntos del intervalo I y obtener una aproximación, y_i, del valor exacto \varphi(x_i), para cada i=0,1, \ldots, N. Supondremos que los puntos son equidistantes, es decir, 

4.6.1. Método de Euler (o de la tangente):

Conocido y_0, el Método de Euler aproxima \varphi(x_1) mediante y_1 \equiv el valor en x_1 de la recta tangente a la curva y=\varphi(x) en el punto (x_0, y_0), cuya pendiente es \varphi'(x_0) = f(x_0,y_0). Esta recta tangente tiene por ecuación

de donde

\varphi(x_1) \simeq y_1 = y_0 + f(x_0,y_0) (x_1 - x_0) = y_0 + h \cdot f(x_0,y_0) .

El proceso se repite, obteniendo:

\varphi(x_2) \simeq y_2 = y_1+ h \cdot f(x_1,y_1) ,

\varphi(x_3) \simeq y_3 = y_2+ h \cdot f(x_2,y_2) 

...

\varphi(x_n) \simeq y_n= y_{n-1}+ h \cdot f(x_{n-1},y_{n-1}) .

Así, el esquema iterativo (o algoritmo) del Método de Euler es:

  (Euler) \quad   \left\{ \begin{array}{l} \mbox{Conocido } y_0,  \\ y_i = y_{i-1}+ h \cdot f(x_{i-1},y_{i-1}), \mbox{ para } i=1, \ldots, N. \end{array} \right.  

Para resolver con Sage el (PVI) mediante el Método de Euler, se utiliza la orden eulers_method, cuya sintaxis es:

O bien, con decimales:

Observamos que Sage devuelve tres columnas, con los valores respectivos de x_i, y_i y h \cdot f(x_i,y_i), para i=0,1,...,N (siendo N=2 en este caso).

o bien con números racionales

Aunque Sage tiene implementado el Método de Euler, también puede programarse. Presentamos a continuación una propuesta de programa y lo aplicamos al ejemplo anterior.

En la primera columna aparecen los valores i=0,1,\ldots,N (que marcan el número de la iteración en la que estamos); en la segunda, los valores x_i (para i=0,1,\ldots,N); y, en la tercera, las aproximaciones y_i (para i=0,1,\ldots,N) de los valores \varphi(x_i) obtenidos por el Método de Euler que acabamos de programar.

También podemos mostrar las salidas con decimales:

4.6.2. Método de Heun (o Método de Euler modificado) 

Conocido y_0, el Método de Heun propone aproximar \varphi(x_1) utilizando, en lugar de la recta tangente a la curva y = \varphi (x) en el punto (x_0,y_0), la recta que pasa por el punto (x_0,y_0)  que tiene por pendiente la media aritmética de las pendientes de las rectas tangentes a y = \varphi (x) en los puntos (x_0, y_0) y (x_1, \widehat{y}_1) , donde  \widehat{y}_1 = y_0 + h \cdot f(x_0, y_0) (el que se obtiene a partir de y_0 mediante el Método de Euler). Es decir, mediante el Método de Heun, se aproxima \varphi(x_1) por y_1 \equiv el valor en x_1 de la recta que pasa por el punto (x_0,y_0) y tiene pendiente 

El proceso se repite para i=1, \ldots, N, obteniendo el siguiente esquema iterativo (o algoritmo) del Método de Heun:

  (Heun) \quad   \left\{ \begin{array}{l} \mbox{Conocido } y_0,  \\  \displaystyle y_i = y_{i-1}+ h \cdot \frac{ f(x_{i-1},y_{i-1}) + f(x_i,  \widehat{y}_i) }{2}, \quad \mbox{para } i=1, \ldots, N, \quad  \mbox{donde} \quad  \widehat{y}_i = y_{i-1}+ h \cdot f(x_{i-1},y_{i-1}). \end{array} \right.  

Una propuesta de programación con Sage del Método de Heun es la siguiente:

Lo aplicamos al método anterior, obteniendo:

O bien con cifras decimales:

4.6.3. Errores de aproximación mediante los Métodos de Euler y de Heun

Puede demostrarse:

• Bajo ciertas condiciones (que f(x,y) posea derivadas parciales continuas), el error local cometido por el Método de Euler en la n-ésima iteración es de orden O(n h^2) (para cada n=1,\ldots,N). Más precisamente,

  donde

• Bajo ciertas condiciones (que f(x,y) posea derivadas parciales segundas continuas), el error local cometido por el Método de Heun en la n-ésima iteración es de orden O(n h^3)   (para cada n=1,\ldots,N). Más precisamente,

Puesto que el paso h>0 está destinado a ser pequeño, el Método de Heun supone una mejora con respecto al Método de Euler.

Por otro lado, si se conoce la solución exacta \varphi(x) del (PVI), se pueden calcular los errores locales exactos

Podemos modificar ligeramente los anteriores programas realizados con Sage para el cálculo de las aproximaciones mediante los Métodos de Euler y de Heun e incluir en la salida una cuarta columna que contenga los errores locales cometidos en cada iteración.

A continuación damos una propuesta de programación con Sage del Método de Euler mostrando las aproximaciones y_i (i=0,1,\ldots,N) y los errores locales correspondientes, y lo aplicamos al Ejemplo 11.

Análogamente, presentamos una propuesta de programación con Sage del Método de Heun mostrando las aproximaciones y_i (i=0,1,\ldots,N) y los errores locales correspondientes, y lo aplicamos al Ejemplo 11.

Finalmente, presentamos dos animaciones. La primera muestra gráficamente la solución exacta, \varphi(x) y la poligonal que une los puntos (x_i,y_i) (para i=0,1,\ldots,N), donde cada y_i es la aproximación de \varphi(x_i) que proporciona el Método de Euler. Nótese que en la ventana interactiva hay que introducir el segundo miembro de la EDO, f(x,y), y la solución exacta del (PVI), \varphi(x).