Orbita plana de una luna entorno a un planeta
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Orbita plana de una luna entorno a un planeta. Grupo 13-B |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2013-14 |
| Autores | Noemí Ortiz, Alicia Chacón, Miguel Sánchez, Mónica Gómez, Cristina Jiménez |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
El trabajo realizado consiste en el estudio de la órbita lunar alrededor de un planeta.
Para ello disponemos de dos ecuaciones diferenciales de segundo grado, que describen dicha trayectoria, y que deberemos convertir en un sistema de cuatro ecuaciones diferenciales no lineales de primer orden, para poder aplicar los métodos de cálculo numérico correspondientes.
\[\left\{\begin{matrix}\ x=-G \frac{mx}{(x^{2}+y^{2})^{3/2}}\ , & t\in [0,T]\\ y=-G \frac{my}{(x^{2}+y^{2})^{3/2}}\ , & \end{matrix}\right.\]
1 Reducción al sistema
Al disponer de una parametrización de dos variables, x(t) e y(t), para transformarlo en el sistema anteriormente citado, hemos tomado dichas variables y sus derivadas, asignando a cada una, una nueva variable.
- [math] x(t)=y_1(t) [/math]
- [math] x'(t)=y_2(t) [/math]
- [math] y(t)=y_3(t) [/math]
- [math] y'(t)=y_4(t) [/math]
De manera que nuestras ecuaciones iniciales quedarían definidas de la siguiente forma:
- [math] y'_2= x''= -G \frac{mx}{(x^{2}+y^{2})^{3/2}}[/math]
- [math] y'_4= y''= -G \frac{my}{(x^{2}+y^{2})^{3/2}}[/math]
Y el sistema final con el que trabajaríamos sería el siguiente: \[\left\{\begin{matrix}\ y'_1(t)=y_2(t) \\ y'_2(t)= -G \frac{my_1}{(y_1^{2}+y_2^{2})^{3/2}}\\ y'_3(t)= y_4(t)\\ y'_4(t)= -G \frac{my_3}{(y_1^{2}+y_3^{2})^{3/2}} & \end{matrix}\right.\]
A continuación explicamos los dos métodos de resolución empleados para este sistema.
2 Método de Euler modificado
Se llama método de Euler al método numérico consistente en ir incrementando paso a paso la variable independiente y hallando su imagen con la derivada correspondiente. Por otra parte, el método de Euler modificado se basa en la misma idea pero hace un refinamiento en la aproximación, tomando un promedio entre ciertas pendientes. Este método es equivalente a un Runge-Kutta de orden 2.
\[y_{n+1}= y_{n}+h\cdot \frac{[(f(x_{n},y_{n}+f(x_{n+1},y_{n+1})]}{2}\] Donde \[y_{n+1}= y_{n}+h\cdot f(x_{n},y_{n})\] Dicho método en código Matlab es el siguiente:
%Médoto Euler modificado
clear all
t0=0;
tf=100;
%longitud de paso
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
% definimos a b c d
a=1; b=0; c=0; d=1;
y0=[a ;b;c;d];
y(:,1)=y0;
yy=y0;
Gm=1;
%Bucle
for n=1:N
r=(yy(1)^2+yy(3)^2)^(3/2);
K1=[yy(2); (-Gm*yy(1)/r);yy(4); (-Gm*yy(3)/r)];
r1=[((yy(1)+h*K1(1))^2+(yy(3)+h*K1(3))^2)^(3/2)];
K2=[yy(2)+h*K1(2); -(Gm*(yy(1)+h*K1(1)))/r1; yy(4)+h*K1(4);-(Gm*(yy(3)+h*K1(3)))/r1];
yy=yy+(h/2)*(K1+K2);
y(:,n+1)=yy;
end
plot(y(1,:),y(3,:))
Además hemos ido variando la longitud de paso para demostrar que la precisión del método es menor cuanto mayor es la longitud de paso.
A continuación se muestran las gráficas correspondientes a dichas variaciones.
3 Método Runge-Kutta
Los métodos de Runge-Kutta son un conjunto de métodos iterativos (implícitos y explícitos) para la aproximación de soluciones de ecuaciones diferenciales ordinarias. En nuestro caso utilizaremos el Runge-Kutta de orden 4 dado por la siguiente ecuación:
- [math]y_{i+1} = y_i + {1 \over 6}h\left ( k_1 + 2k_2 + 2k_3 + k_4 \right ) [/math]
Donde
- [math] \begin{cases} k_1 & = f \left( x_i, y_i \right) \\ k_2 & = f \left( x_i + {1 \over 2}h, y_i + {1 \over 2}k_1 h \right) \\ k_3 & = f \left( x_i + {1 \over 2}h, y_i + {1 \over 2}k_2 h \right) \\ k_4 & = f \left( x_i + h, y_i + k_3h \right) \\ \end{cases} [/math]
A continuación mostramos el código Matlab utilizado para la descripción de la trayectoria utilizando este método:
%metodo Runge-Kutta de cuarto orden
clear all
t0=0;
tf=2000;
%longitud de paso
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%Definimos a, b, c y d
a=1; b=0; c=0; d=1;
yo=[a;b;c;d]
y(:,1)=yo;
yy=yo;
Gm=1.1;
%Empezamos el bucle
for n=1:N
r=(yy(1)^2+yy(3)^2)^(3/2);
K1=[yy(2); (-Gm*yy(1)/r);yy(4); (-Gm*yy(3)/r)];
r1=[((yy(1)+h*K1(1))^2+(yy(3)+h*K1(3))^2)^(3/2)];
K2=[yy(2)+(h/2)*K1(2); -(Gm*(yy(1)+(h/2)*K1(1)))/r1; yy(4)+(h/2)*K1(4);-(Gm*(yy(3)+(h/2)*K1(3)))/r1];
K3=[yy(2)+(h/2)*K2(2); -(Gm*(yy(1)+(h/2)*K2(1)))/r1; yy(4)+(h/2)*K2(4);-(Gm*(yy(3)+(h/2)*K2(3)))/r1];
K4=[yy(2)+(h)*K3(2); -(Gm*(yy(1)+(h)*K3(1)))/r1; yy(4)+(h)*K3(4);-(Gm*(yy(3)+(h)*K3(3)))/r1];
yy=yy+(h/6)*(K1+2*(K2)+2*K3+K4);
y(:,n+1)=yy;
end
plot(y(1,:),y(3,:))
Y al igual que en el método anterior, hemos obtenido diferentes gráficas en función de los diferentes valores que hemos asignado a la longitud de paso y a los extremos del intervalo:
