Metodo Euler Y Runge-Kutta 4th Order Para EDO De Orden Superior

El propósito del siguiente proyecto es poder guiar al usuario a aprender algunos conceptos básicos de como resolver ecuaciones diferenciales de orden superior usando métodos numéricos: Método de Euler y Runge-Kutta 4to Orden. Esta guía está diseñada para que el usuario, incluso sin experiencia en la utilizacion de metodos numericos, pueda de aquí en adelante llevar a cabo este tipo de operaciones. Para ello, a continuacion podran ver los pasos a seguir para llevar a cabo cada ejemplo. Vale mencionar que los métodos numéricos se ejecutan usando el programa MatLab el cual tiene precio, sin embargo, el programa Octave es gratis y corre el mismo tipo de archivo. La diferencia que si hay que tomar en cuenta es que los comentarios del programa, en matlab, se escriben con el uso de "%" a principios de la fila. En contraste, en Octave se utiliza el símbolo "$" para referirse a un comentario

• Matlab/Octave

Matlab es un programa cuyo rango de aplicación es inmenso por lo cual se le recomienda a los estudiantes de ingeniería. Con este programa se puede graficar, resolver e incluso diseñar sistemas de fenómenos físicos para así analizar un prototipo. Sin embargo, en este caso se utilizó para el resolver Ecuaciones diferenciales de orden superior. Para ello se espera que el usuario conozca los siguiente temas: ciclos "for", vectores y graficación. Por ello a continuación se darán una introducción a lo que lleva a cabo cada comando.

Ciclos "for"

Los ciclos for se utilizan para ejecutar una acción de forma continua. Para ello se debe establecer condiciones iniciales y condiciones finales. Ya que al definirse se establece que "para" cada valor de cierta variable que cumpla con los parámetros, el programa ejecutará ciertas acciones.

Vectores

Los vectores sirven como contenedores de almacenamiento. A partir de ello se le puede indicar al código que se desea llamar al elemento (x) del vector. Vale mencionar que en Matlab, a diferencia de Python, el primer elemento del vector tiene asignado el (1) y a partir de ello se sigue creciendo (en Python se inicia con 0).

Graficación

La graficación de los valores presenta posibilidades de graficar datos en 2 dimensiones al igual que en 3 dimensiones. Este tipo de acciones se logra al utilizar el comando plot() y dentro del paréntesis las entradas se separan por una coma(,). Vale mencionarse que en el paréntesis el orden de lo que se ha de llamar va ser el siguiente:

1. valor que se va a graficar en el eje x
2. valor que se va a graficar en el eje y
3. valor que se va a graficar en el eje z (OPCIONAL)
4. característica de la gráfica ya sea apariencia o color (OPCIONAL)

Ej. plot(x,y,z, 'r*')

Para este trabajo, en ambos métodos solo se graficó los valores en 2D por lo que no hay necesidad de incluir la entrada número 3.

Se llama método de Euler al método numérico consistente en ir incrementando paso a paso la variable independiente y hallando la siguiente con la derivada. Usando la fórmula de arriba se predice un nuevo valor de Y por medio de la pendiente (igual a la primera derivada en el valor original de X). Tomando esto en cuenta vale mencionar que el Método de Euler solo funciona con ecuaciones diferenciales de primer orden, Sin embargo a partir de sustitución se puede convertir una ecuación diferencial de orden superior a una EDO de primer orden.

El ejemplo que se utilizó para este método es una sistema de masa y resorte sin amortiguamiento. Se reconoce que a partir de la ley de Newton, Fuerza es igual a la masa por la aceleración. Sin embargo, en un sistema de resorte, se sabe que la Fuerza ejercida por el resorte es equivalente a su desplazamiento y proporcional a la constante negativa de este. Teniendo ya esto se establece la relación:

ma = -kx

A partir de esto se puede empezar a hacer sustituciones dentro de la misma ecuación. Se reconoce que la derivada respecto al tiempo de la posición o la pendiente (dx/dt) es la velocidad en un punto específico. Ademas de ello, la segunda derivada de la posición respecto al tiempo (d^2x/dt^2) es equivalente a la aceleración (a). Usando esto se puede llegar a formar la siguiente fórmula:

m*d^2x/dt^2 = -kx

Sin embargo, Euler no es capaz de operar si se utiliza derivadas de segundo orden o mayor. Sin embargo aquí es donde se puede hacer la sustitución. Ya que se conoce que la velocidad es la primera derivada de la posición y la aceleración es la segunda, se puede decir que la aceleración es la primera derivada de la velocidad respecto al tiempo (dv/dt). Tomando esto en cuenta y luego despejando la derivada, se puede llegar a la siguiente relación:

m*dv/dt = -kx

dv/dt = -kx/m

A partir de intuición se puede reconocer que el comportamiento de éste será oscilatorio. En nuestro modelo, no hay algún tipo de amortiguamiento por lo que la energía se conserva y su posición inicial en t=0 es y_0 = 1. Ademas de ello, en nuestro programa se consideró que el coeficiente del resorte y la masa fuesen equivalentes. Por ello la ecuación real del sistema que modelaría su posición en el eje y sobre el tiempo sería la siguiente fórmula:

y = y_0*cos (t*k/m)

y = cos(t)

Como se puede observar en el código, dentro del ciclo "for" se encuentra la ecuación que genera el comportamiento real, la cual está denominada como "f_1." Adicionalmente, la formula de Euler (La cual se puede apreciar en la imagen superior) se replico consecuentemente dentro del ciclo. Esto asegura que el programa se adaptara a las condiciones anteriormente especificadas o predecidas, lo cual genera la gráfica de Euler.

El método de Runge-Kutta no es solo un único método, sino una familia de métodos iterativos, tanto implícito como explicito, para aproximar las soluciones de ecuaciones diferenciales ordinarias. El más usado de la familia es Runge-Kutta de cuarto orden, a menudo es referido como RK4. Este método se define por las ecuaciones que se pueden apreciar en la parte superior. Para nuestro modelo se buscó la graficación de la posición de una partícula al igual que su posición respecto al tiempo.

De la misma forma que el método de Euler, RK4 solo puede funcionar con derivadas de primer orden, para lograr esto se ha de sustituir cualquier derivada de orden superior para así formar un sistema de primeras derivadas. para empezar se define que la aceleración del sistema es la derivada de la velocidad, por lo que se puede identificar en el código como "fv." De igual forma, la derivada de la posición es su velocidad, por lo que la velocidad se puede definir inicialmente como "fx." Nuestro periodo de tiempo se definirá como "t" el cual será desde un punto inicial "t0" a un punto final "tmax" en intervalos de "h."

Dentro del ciclo, como se puede ver en la imagen, se empieza a hacer la definición y cálculo de los distintos valores de velocidad y posición. Esto se llevará a cabo durante un periodo de tiempo definido. Por ultimo, vale mencionar que después del ciclo for se puede ver lo que sería la ejecución de graficación del sistema. Aquí es donde se define qué se graficaron las distintas figuras, se definirá el nombre de cada una, el titulo de la grafica, y entre que intervalo se enseñara los resultados.

Para poder alterar nuestro sistema para futuras ocasiones solo se ha de hacerle cambios a los valores de posicion.y velocidad que el usuario desee. Esto se debe a que el método está definido dentro del ciclo. Ahora si se desea evitar el paso de graficar o incluso imprimir los valores de velocidad y posición conforme pasa el tiempo, pues se puede hacer los comandos comentarios al ponerle un "%" previo a la declaración en la misma fila.