II. Solución numérica de la ecuación del péndulo simple.

No había hecho problemas como este, si que había resuelto ecuaciones diferenciales pero no por métodos numéricos. Asi que todo lo que he hecho lo encontré en esta página:

<a>http://galia.fc.uaslp.mx/~jvallejo/EDO.pdf</a>

Ya que la ayuda de Máxima era insuficiente. Y bueno, yo te voy a enseñar con un programa que a lo mejor es menos amigable que otros (de hecho lo es), pero es con el único que me manejo un poco.

En la página 82 y anyeriores de ese manual tenemos lo que nos interesa . Usaremos el método de Runge Kutta de orden 4 que hay en Máxima. En el ejemplo, se resuelve una ecuación y se compara con la respuesta exacta. Aquí nos dicen que comparemos con la respuesta de un sistema aproximado. Luego en vez de la gráfica de la ecuación exacta, usaremos la de este sistema aproximado.

Asignando el valor a las constantes el sistema de ecuaciones es:

$$\begin{align}&\frac{dx_1}{dt}=x_2\\ &\\ &\frac{dx_2}{dt}=-49 sen(x_1)-\frac{x_2}2\\ &\\ &\text{Y el sistema aproximado que dicen será:}\\ &\\ &\frac{d^2x_1}{dt^2}=-49x_1 -\frac 12 \frac{dx_1}{dt}\\ &\end{align}$$

Te escribo lás órdenes y si no entiendes algo puedes consultarlo en el manual.

Primero estas para cargar los paquetes necesarios

load(dynamics);

load(draw);

Ahora estas otras que calculan la solución exacta del sistema aproximado con el que nos piden comparar

ec: 'diff(x1,t,2)= -49*x1 - (1/2)*'diff(x1,t);

solge: ode2(ec,x1,t);

ic2(solge, t=0, x1=%pi/3, 'diff(x1,t)=0);

Saldrá una expresión bastante complicada

x1=%e^(-t/4)*((%pi*sin((3*sqrt(87)*t)/4))/(9*sqrt(87))+(%pi*cos((3*sqrt(87)*t)/4))/3)

Que no es cuestión de que tengamos que escribirla a mano, se puede seleccionar y copiar para usar en otros sitios pegándola.

Ahora se calcula la tabla de valores del método Runge-Kutta. Por lo que sea no funciona si se escribe %pi/3, por eso he puesto el valor numérico de eso

tabla: rk([x2,-49*sin(x1)-0.5*x2],[x1,x2],[1.047197551,0],[t,0,5,0.05]);

Se crea una segunda tabla que solo tiene las coordenadas x1, x2 sin la t

tablab: makelist([tabla[1],tabla[2]],i,1,length(tabla));

Y ahora graficamos la solución exacta del sistema aproximado que nos dicen y la tabla Runge-Kutta de la resolución numérica del sistema. En la siguiente orden es donde pegamos en un sitio la solución que habíamos obtenido arriba para el sistema aproximado.

wxdraw2d(key="sol. exacta aproximada", explicit(%e^(-t/4)*((3.141592654 * sin((3*sqrt(87)*t)/4))/(9*sqrt(87))+(3.141592654 * cos((3*sqrt(87)*t)/4))/3), t, 0, 5), points_joined = true,color="red", key="Runge-Kutta con h=0.05", points(tablab));

Y el resultado es este:

He hecho la gráfica con un paso más pequeño h=0.001 y la solución Runge-Kutta es más o menos la misma. Yo no me fio mucho de que el sistema que ne nos han dado para comparar sea tan aproximado al real. Pero lo que pasa es que el sistema real no me lo resuelve ni Máxima noi Wolframalpha, luego no puedo asegurar nada.

Me ha costado mucho. Estuve intentando en vano que la solución se pudiese transformar en función para no tener que hacer ese copiar y pegar, para ello miré en muchos sitios pero no encontré respuesta. Te mando esto de momento. Si me mandas otra pregunta para lo que queda del ejercicio estaría encantado.

Un saludo.

Pregunté personalmente al responsable de la traducción al español de la ayuda de Maxima y he obtenido la respuesta a lo de no tener que copiar la expresión.

Las órdenes serán estas:

load(dynamics);

load(draw);

ec: 'diff(x1,t,2) = -49*x1 - (1/2)*'diff(x1,t);
solge: ode2(ec,x1,t);
solpar:rhs( ic2(solge, t=0, x1=%pi/3, 'diff(x1,t)=0));

tabla: rk([x2,-49*sin(x1)-0.5*x2],[x1,x2],[1.047197551,0],[t,0,5,0.05]);

tablab: makelist([tabla[i][1],tabla[i][2]],i,1,length(tabla));

wxdraw2d(key="sol. exacta aproximada", explicit(solpar, t, 0, 5), points_joined = true,color="red", key="Runge-Kutta con h=0.05", points(tablab));

Se ha modificado una línea por esta otra:

solpar:rhs( ic2(solge, t=0, x1=%pi/3, 'diff(x1,t)=0));

Que es la que hace que solpar sea una expresión que representa a la solución como función de t

Y se ha corregido esta otra

tablab: makelist([tabla[i][1],tabla[i][2]],i,1,length(tabla));

donde habían desaparecido unos [i] de manera misteriosa, no se si me los olvidé o los suprimió la página.

Y ahora si que debe funcionar, yo he copiado esas órdenes de auí a Maxima y ha funcionado.

Que asco, han hecho más estrecha y menos alta la ventana de escritura, no aprovecha ni la mitad del ancho del navegador, pero los facebook, twitter y otras cosillas no pueden faltar.

Me ha llegado un correo de Todo Expertos como que me has hecho una nueva pregunta, pero no puedo llegar a ella, si sigo el enlace del correo me lleva a un sitio donde no sale nada. Ha habido cambios muy radicales en la página, prueba a hacer la pregunta de nuevo, mándala a la página en vez a mí, no sé, intenta algo porque de la forma que lo has hecho no aparece la pregunta por ningún lado.

No perdona! No era un mensaje nuevo, era uno viejo pero me ha parecido que era recién recibido. Y ya no se puede acceder a cosas que existían hace dos días con el mismo enlace, lo han cambiado todo.