ER1.7. Solución numérica de ecuaciones diferenciales

Introducción (planteamiento)

Dadas las siguientes ecuaciones diferenciales no lineales, resolverlas por el método de Euler con un paso (tiempo de muestreo) suficientemente pequeño, comparar con los resultados utilizando la función ode45 de MATLAB (basada en una fórmula explícita de Runge-Kutta (4,5), el método de Dormand-Prince) y analizar los resultados:

1. $\ddot{y}+3\dot{y}+2y=1, y(0) =1, \dot{y}(0) =1$
2. $\frac{dy}{dt}=t^2+y^2, y(0)=1$

Método (plan de solución)

Pasos: (i) si la ecuación diferencial es de orden mayor que dos, se debe convertir a una ecuación de estado; (ii) implementar en MATLAB el método de Euler con diferentes pasos; (iii) utilizar la función ode45 de MATLAB; (iv) interpretar los resultados y comparar el primer caso con los resultados del ejemplo 1.6.

Resultados (solución)

1. $\ddot{y}+3\dot{y}+2y=1, y(0) =1, \dot{y}(0) =1$

Lo primero es pasar la ecuación diferencial a la forma en variables de fase:

$x_1=y, x_2=\dot{y}$ 

$\begin{cases} \dot{x}_1=x_2\\ \dot{x}_2=-2x_1-3x_2+1\\\end{cases}$ 

$x_1(0)=1, x_2(0)=1$

Aplicando el método de Euler en diferencias hacia delante se obtiene:

$\begin{cases} x_1(k+1)=x_1(k)+T_sx_2(k) \\ x_2(k+1)=x_2(k)+T_s\left[ -2x_1(k)-3x_2(k)+1 \right]\\ \end{cases}$

El código de MATLAB de este enlace implementa las ecuaciones anteriores con $T_s=1$ y compara las soluciones numéricas con la solución analítica y la solución numérica por el método ode45. En el ejemplo resuelto 1.6 se ilustra el cambio del período de muestreo en este ejemplo.

Resultado:

Los resultados numéricos coinciden con la solución analítica y, adicionalmente, se obtiene la solución $x_2=\dot{y}$ (se invita al lector a revisar si la segunda curva realmente es la derivada de la primera).

2. $\frac{dy}{dt}=t^2+y^2, y(0)=1$

El método de Euler da la siguiente ecuación en diferencias (la cual se puede resolver iterativamente), donde es necesario hacer el siguiente cambio explícito en $t=kT_s$, pero no en $y(t)$:

$y(k+1)=y(k)+T_s\left[( kT_s) ^2+y^2(k) \right] , y(0)=1$

La ecuación anterior para diferentes pasos y el método ode45 se implementan con el código de MATLAB de este enlace.

Resultado: 

Discusión y verificación

Los ejemplos anteriores muestran que el método de Euler entrega resultados aceptables si el período de muestreo es lo suficientemente pequeño. Sin embargo, el método de Runge-Kutta implementado en la función ode45 de MATLAB da buenos resultados con un paso mayor, por lo que debe ser el método por defecto. En todos los casos, se debe convertir la ecuación diferencial a una ecuación de estado en variables de fase.

El segundo ejercicio muestra que aunque se cumplen las condiciones del teorema de existencia y unicidad (el lector lo puede comprobar), no se sabe fácilmente para que intervalo de la variable independiente la solución existe. Si se toma un tiempo final mayor que 0.96 la solución deja de existir (se invita al lector a verificar a partir del código dado).