Resolviendo ecuaciones diferenciales ordinarias de primer orden con el método de Euler

 El método de Euler no es precisamente el mejor; pero para si es lo más recomendable para iniciar con la programación de este tipo de métodos. Hay otros métodos con mejor precisión y robustos; pero ya los iremos abordando partiendo de lo que aprendamos de la implementación del método de Euler.

Primero, hay unas cosas que debemos saber de este método:

• Solo funciona para ecuaciones diferenciales ordinarias de primer orden
• No importa si la ecuación es de coeficientes constantes o de coeficientes variables
• Si la ecuación de un orden mayor, debes usar tantos cambios de variables como sean necesarios

No vamos a explicar las bases del método de Euler; pero sí les puedo recomendar algunas fuentes de fácil acceso:

• Métodos Numéricos para Ingenieros. S. C. Chapra y R. P. Canale. Yo tengo las ediciones 3 y 7, y la explicación es bastante didáctica. Además tiene ejemplos que puedes poner a prueba en MAPLE
• Ecuaciones Diferenciales con Aplicaciones de Modelado. D. G. Zill. Este libro no es complicado y tiene un apartado de métodos numéricos para ecuaciones diferenciales. También es fácil de entender

En resumen, el método de Euler establece que puedes resolver numéricamente una ecuación diferencial como la siguiente:

$\frac{dy}{dx}=f(x,y)$

usando la siguiente fórmula:

$y_{i+1}=y_i+f\left( x_i, y_i \right) h$

en donde $h$ es el tamaño de paso y la $i$ es un contador numérico entero que depende de $h$ y del tamaño del intervalo en donde se desee resolver la ecuación diferencial. La fórmula anterior es una versión discretizada de la ecuación original.

La cuestión importante aquí es que si nosotros conocemos el valor de $y$ en algún valor de $x$, entonces podríamos calcular otro valor de $y$ que corresponda a otro valor de $x$. En este punto el rol que tiene la $i$ es importante porque define o etiqueta si un valor de $x$ y $y$ es anterior o posterior en el cálculo. Dicho de otro modo, los $y_{i+1}$ se calculan a partir de los $y_i$ y de los $x_i$.

Un comando en MAPLE

En realidad no es un comando, sino una opción de un comando. MAPLE tiene una opción para resolver ecuaciones diferenciales ordinarias con muchos métodos, y el de Euler es uno de estos. Consideremos entonces la siguiente ecuación diferencial:

$\frac{dy}{dx}=-2x^3+12x^2-20x+8.5$

sujeta a la siguiente condición inicial:

$y=1$ en $x=0$

Este ejemplo lo tomé del libro de Chapra y Canale solo para mostrar que llegamos al mismo resultado. La implementación en MAPLE es como sigue:

> with(Student[NumericalAnalysis]):

> ode1:= diff(y(x),x)=-2*x**3+12*x**2-20*x+8.5;

> InitialValueProblem(ode1, y(0) = 1, x = 4, method = euler, output = information, numsteps = 8, digits = 6);

En este caso, se tomó $h=0.5$. La tercera columna de los resultados anteriores muestra lo que uno obtendría si hiciera los cálculos a mano. Lo mejor de todo es que también proporciona la solución que MAPLE obtiene numéricamente (con su mejor método). Por supuesto que también en este caso es fácil obtener la solución exacta de la ecuación diferencial problema integrando una sola vez y comparar con la tabla anterior.

La cuarta columna muestra porque el método de Euler no es tan usado. El error aumenta conforme nos alejamos del punto inicial. Sin embargo, si disminuimos el tamaño de paso los resultados mejoran. Por ejemplo, entre $x=0$ y $x=2$, y con un tamaño de paso $h=0.25$ se obtiene:

Si comparamos los resultados para $x=2$ vemos que el resultado con el método de Euler mejora.

Te recomiendo que busques en la ayuda InitialValueProblem para que tengas una idea más amplia de como usar este comando.

Implementando el método de Euler en MAPLE

Programar el método de Euler en MAPLE no es complicado. Chapra y Canale muestran dos algoritmos posibles. Aquí abordaremos la forma más sencilla implementando una estructura for. El código sería como sigue:

> y[1]:= 1:

  x[1]:= 0:# Condición inicial

> hstp:= 0.5:# Tamaño de paso

> linf:= x[1]:# Inicio del intervalo

  lsup:= 4.0:# Fin del intervalo

> cnt:= floor( (lsup - linf)/hstp ):# Contador 

> fii:= - 2*x[i]**3 + 12*x[i]**2 - 20*x[i] + 8.5:# Inhomogeneidad de ODE

> # Cálculos del método de Euler

  print("x", "y");

  print(x[1], y[1]);

  for i from 1 to cnt do

    y[i+1]:= y[i] + fii*hstp:

    x[i+1]:= x[i] + hstp:  

    print(x[i+1], y[i+1]);

  end do:

  unassign('i');

La implementación es muy simple y la presentación de resultados quizá no es muy estética; pero puede comprobarse que coincide con lo que se obtiene con el comando de MAPLE de la sección anterior. Para el contador utilicé el comando floor para que se usará siempre un número entero; aunque en realidad MAPLE no tendría problema si le proporcionamos un contador no entero en la estructura for.

Debe notarse como intenté que el programa tuviera cierta adaptabilidad a otros problemas o requerimientos. Por ejemplo, no es dificil:

• cambiar el tamaño de paso,
• cambiar el tamaño del intervalo,
• si cambiamos de ecuación solo debes actualizar fii.

Si tienes alguna duda escríbela en los comentarios y trataré de ayudarte.

Ildebrando