Diferencia entre revisiones de «Órbita lunar»

De MateWiki
Saltar a: navegación, buscar
Línea 72: Línea 72:
 
plot(x1,y1,'r')
 
plot(x1,y1,'r')
  
%para h=1
 
 
xx1=zeros(N2+1);
 
xx1(1)=1;
 
xx2=zeros(N2+1);
 
xx2(1)=0;
 
yy1=zeros(N2+1);
 
yy1(1)=0;
 
yy2=zeros(N2+1);
 
yy2(1)=1;
 
Z1=zeros(4,N2+1);
 
 
for n=1:N2
 
    ZZ1=[xx1(n);xx2(n);yy1(n);yy2(n)];
 
    kk1=[ZZ1(2);-((ZZ1(1))/(((ZZ1(1))^2)+(ZZ1(3)^2))^(3/2));ZZ1(4);-((ZZ1(3))/(((ZZ1(1))^2)+...
 
        (ZZ1(3)^2))^(3/2))];
 
    Zp1=ZZ1+h2*kk1;
 
    kk2=[Zp1(2);-((Zp1(1))/(((Zp1(1))^2)+(Zp1(3)^2))^(3/2));Zp1(4);...
 
        -((Zp1(3))/(((Zp1(1))^2)+(Zp1(3)^2))^(3/2))];
 
    Z1(:,n+1)=ZZ1+(h2*((kk1+kk2)/2));
 
    xx1(n+1)=Z1(1,n+1);
 
    xx2(n+1)=Z1(2,n+1);
 
    yy1(n+1)=Z1(3,n+1);
 
    yy2(n+1)=Z1(4,n+1);
 
end
 
 
figure(2)
 
clf
 
plot(xx1,yy1,'g')
 
 
}}
 
}}
  

Revisión del 15:50 28 feb 2014

Trabajo realizado por estudiantes
Título Órbita lunar. Grupo 5-B
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores Francisco Durán Muñoz, Javier Bosch Martínez, Manuel López Martín, Miguel García García, Emilio Valero Muñoz-Rojas, Javier Tordera Garrido
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Rotación es el movimiento de cambio de orientación de un cuerpo o un sistema de referencia de forma que una línea (llamada eje de rotación) o un punto permanece fijo. La posición estaría definida por las coordenadas (x,y), que dependen de una función. En nuestro caso conocemos las aceleraciones:


\[\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 de un sistema

Consiste en reducir una ecuación de orden n, a n ecuaciones de orden 1. Como en este caso tenemos un sistema de 2 ecuaciones de orden 2, nos quedará un nuevo sistema formado por 4 ecuaciones de orden 1. Realizando el siguiente cambio de variable:

[math]x1  =  x \quad  x2 = x'[/math]
[math]y1  =  y \quad  y2 = y'[/math]


Nos quedarán las siguientes fórmulas: \[x1'=x2\] \[y1'=y2\]

\[x2'=-G\frac{m\cdot x1}{(x1^{2}+y1^{2})^{3/2}}\] \[y2'=-G\frac{m\cdot y1}{(x1^{2}+y1^{2})^{3/2}}\]

Tierra.gif

2 Método de Euler modificado

El método de Euler se basa en aproximar el valor de la función a la tangente en cada punto, cuanto más cercanos cojamos los puntos, mayor será la aproximación al resultado real, haciéndolo más exacto .El método de Euler modificado es utilizado para resolver ecuaciones diferenciales no lineales, se basa en la misma idea del método anterior, pero hace un refinamiento en la aproximación, tomando un promedio entre ciertas pendientes.

La fórmula es la siguiente: \[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})\] El programa quedaría:

%trabajo apartado 2
clear all
t0=0;
tN=100;
h1=0.1;
h2=1;
N1=(tN-t0)/h1;
N2=(tN-t0)/h2;
t1=t0:h1:tN;
t2=t0:h2:tN;

%para h=0.1

x1=zeros(N1+1);
x1(1)=1;
x2=zeros(N1+1);
x2(1)=0;
y1=zeros(N1+1);
y1(1)=0;
y2=zeros(N1+1);
y2(1)=1;
Z=zeros(4,N1+1);
for n=1:N1
    ZZ=[x1(n);x2(n);y1(n);y2(n)];
    k1=[ZZ(2);-((ZZ(1))/(((ZZ(1))^2)+(ZZ(3)^2))^(3/2));ZZ(4);-((ZZ(3))/(((ZZ(1))^2)+...
        (ZZ(3)^2))^(3/2))];
    Zp=ZZ+h1*k1;
    k2=[Zp(2);-((Zp(1))/(((Zp(1))^2)+(Zp(3)^2))^(3/2));Zp(4);-((Zp(3))/(((Zp(1))^2)+(Zp(3)^2))^(3/2))];
    Z(:,n+1)=ZZ+(h1*((k1+k2)/2));
    x1(n+1)=Z(1,n+1);
    x2(n+1)=Z(2,n+1);
    y1(n+1)=Z(3,n+1);
    y2(n+1)=Z(4,n+1);
end
figure(1)
clf
plot(x1,y1,'r')


3 Método Runge-Kutta

Los métodos de Runge-Kutta son un conjuntos de métodos iterativos (implícitos y explícitos) para la aproximación de soluciones de ecuaciones diferenciales ordinarias, concretamente, del problema de valor inicial. Sabiendo que el método de Runge-Kutta consiste en: [math]y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)[/math] En nuestro caso, el programa queda de la siguiente manera:

%apartado3 runge kuta
clear all
t0=0;
t1=100;
a1=1;a2=1.3;
b=0;c=0;
d1=1;d2=1.1;d3=1.4;
h=1;
N1=(t1-t0)/h;
x1=zeros(N1+1);
x2=zeros(N1+1);
y1=zeros(N1+1);
y2=zeros(N1+1);
Z1=zeros(4,N1+1);
for n=1:N1
    x1(1)=a1;x2(1)=b;y1(1)=c;y2(1)=d1;
    G=1;
    ZZ=[x1(n);x2(n);y1(n);y2(n)];
    k1=[ZZ(2);-((G*ZZ(1))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2));ZZ(4);...
        -((G*ZZ(3))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2))];
    Zp=ZZ+((h/2)*k1);
    k2=[Zp(2);-((G*Zp(1))/(((Zp(1))^2)+((Zp(3))^2))^(3/2));Zp(4);...
        -((G*Zp(3))/(((Zp(1))^2)+((Zp(3))^2))^(3/2))];
    Zx=ZZ+((h/2)*k2);
    k3=[Zx(2);-((G*Zx(1))/(((Zx(1))^2)+((Zx(3))^2))^(3/2));Zx(4);...
        -((G*Zx(3))/(((Zx(1))^2)+((Zx(3))^2))^(3/2))];
    Zw=ZZ+(h*k3);
    k4=[Zw(2);-((G*Zw(1))/(((Zw(1))^2)+((Zw(3))^2))^(3/2));Zw(4)...
        ;-((G*Zw(3))/(((Zw(1))^2)+((Zw(3))^2))^(3/2))];
    Z1(:,n+1)=ZZ+((h/6)*(k1+(2*k2)+(2*k3)+k4));
    x1(n+1)=Z1(1,n+1);
    x2(n+1)=Z1(2,n+1);
    y1(n+1)=Z1(3,n+1);
    y2(n+1)=Z1(4,n+1);
end

figure(1)
clf
plot(Z1(1,:),Z1(3,:),'r')

x1=zeros(N1+1);
x2=zeros(N1+1);
y1=zeros(N1+1);
y2=zeros(N1+1);
Z1=zeros(4,N1+1);
for n=1:N1
    x1(1)=a2;x2(1)=b;y1(1)=c;y2(1)=d1;
    G=1;
    ZZ=[x1(n);x2(n);y1(n);y2(n)];
    k1=[ZZ(2);-((G*ZZ(1))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2));ZZ(4);...
        -((G*ZZ(3))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2))];
    Zp=ZZ+((h/2)*k1);
    k2=[Zp(2);-((G*Zp(1))/(((Zp(1))^2)+((Zp(3))^2))^(3/2));Zp(4);...
        -((G*Zp(3))/(((Zp(1))^2)+((Zp(3))^2))^(3/2))];
    Zx=ZZ+((h/2)*k2);
    k3=[Zx(2);-((G*Zx(1))/(((Zx(1))^2)+((Zx(3))^2))^(3/2));Zx(4);...
        -((G*Zx(3))/(((Zx(1))^2)+((Zx(3))^2))^(3/2))];
    Zw=ZZ+(h*k3);
    k4=[Zw(2);-((G*Zw(1))/(((Zw(1))^2)+((Zw(3))^2))^(3/2));Zw(4)...
        ;-((G*Zw(3))/(((Zw(1))^2)+((Zw(3))^2))^(3/2))];
    Z1(:,n+1)=ZZ+((h/6)*(k1+(2*k2)+(2*k3)+k4));
    x1(n+1)=Z1(1,n+1);
    x2(n+1)=Z1(2,n+1);
    y1(n+1)=Z1(3,n+1);
    y2(n+1)=Z1(4,n+1);
end

figure(2)
clf
plot(x1,y1,'g')

x1=zeros(N1+1);
x2=zeros(N1+1);
y1=zeros(N1+1);
y2=zeros(N1+1);
Z1=zeros(4,N1+1);
for n=1:N1
    x1(1)=a2;x2(1)=b;y1(1)=c;y2(1)=d2;
    G=1;
    ZZ=[x1(n);x2(n);y1(n);y2(n)];
    k1=[ZZ(2);-((G*ZZ(1))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2));ZZ(4);...
        -((G*ZZ(3))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2))];
    Zp=ZZ+((h/2)*k1);
    k2=[Zp(2);-((G*Zp(1))/(((Zp(1))^2)+((Zp(3))^2))^(3/2));Zp(4);...
        -((G*Zp(3))/(((Zp(1))^2)+((Zp(3))^2))^(3/2))];
    Zx=ZZ+((h/2)*k2);
    k3=[Zx(2);-((G*Zx(1))/(((Zx(1))^2)+((Zx(3))^2))^(3/2));Zx(4);...
        -((G*Zx(3))/(((Zx(1))^2)+((Zx(3))^2))^(3/2))];
    Zw=ZZ+(h*k3);
    k4=[Zw(2);-((G*Zw(1))/(((Zw(1))^2)+((Zw(3))^2))^(3/2));Zw(4)...
        ;-((G*Zw(3))/(((Zw(1))^2)+((Zw(3))^2))^(3/2))];
    Z1(:,n+1)=ZZ+((h/6)*(k1+(2*k2)+(2*k3)+k4));
    x1(n+1)=Z1(1,n+1);
    x2(n+1)=Z1(2,n+1);
    y1(n+1)=Z1(3,n+1);
    y2(n+1)=Z1(4,n+1);
end

figure(3)
clf
plot(x1,y1,'c')

x1=zeros(N1+1);
x2=zeros(N1+1);
y1=zeros(N1+1);
y2=zeros(N1+1);
Z1=zeros(4,N1+1);
for n=1:N1
    x1(1)=a1;x2(1)=b;y1(1)=c;y2(1)=d3;
    G=1;
    ZZ=[x1(n);x2(n);y1(n);y2(n)];
    k1=[ZZ(2);-((G*ZZ(1))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2));ZZ(4);...
        -((G*ZZ(3))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2))];
    Zp=ZZ+((h/2)*k1);
    k2=[Zp(2);-((G*Zp(1))/(((Zp(1))^2)+((Zp(3))^2))^(3/2));Zp(4);...
        -((G*Zp(3))/(((Zp(1))^2)+((Zp(3))^2))^(3/2))];
    Zx=ZZ+((h/2)*k2);
    k3=[Zx(2);-((G*Zx(1))/(((Zx(1))^2)+((Zx(3))^2))^(3/2));Zx(4);...
        -((G*Zx(3))/(((Zx(1))^2)+((Zx(3))^2))^(3/2))];
    Zw=ZZ+(h*k3);
    k4=[Zw(2);-((G*Zw(1))/(((Zw(1))^2)+((Zw(3))^2))^(3/2));Zw(4)...
        ;-((G*Zw(3))/(((Zw(1))^2)+((Zw(3))^2))^(3/2))];
    Z1(:,n+1)=ZZ+((h/6)*(k1+(2*k2)+(2*k3)+k4));
    x1(n+1)=Z1(1,n+1);
    x2(n+1)=Z1(2,n+1);
    y1(n+1)=Z1(3,n+1);
    y2(n+1)=Z1(4,n+1);
end

figure(4)
clf
plot(x1,y1,'k')


x1=zeros(N1+1);
x2=zeros(N1+1);
y1=zeros(N1+1);
y2=zeros(N1+1);
Z1=zeros(4,N1+1);
for n=1:N1
    x1(1)=a1;x2(1)=b;y1(1)=c;y2(1)=d1;
    G=1.1;
    ZZ=[x1(n);x2(n);y1(n);y2(n)];
    k1=[ZZ(2);-((G*ZZ(1))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2));ZZ(4);...
        -((G*ZZ(3))/(((ZZ(1))^2)+((ZZ(3))^2))^(3/2))];
    Zp=ZZ+((h/2)*k1);
    k2=[Zp(2);-((G*Zp(1))/(((Zp(1))^2)+((Zp(3))^2))^(3/2));Zp(4);...
        -((G*Zp(3))/(((Zp(1))^2)+((Zp(3))^2))^(3/2))];
    Zx=ZZ+((h/2)*k2);
    k3=[Zx(2);-((G*Zx(1))/(((Zx(1))^2)+((Zx(3))^2))^(3/2));Zx(4);...
        -((G*Zx(3))/(((Zx(1))^2)+((Zx(3))^2))^(3/2))];
    Zw=ZZ+(h*k3);
    k4=[Zw(2);-((G*Zw(1))/(((Zw(1))^2)+((Zw(3))^2))^(3/2));Zw(4)...
        ;-((G*Zw(3))/(((Zw(1))^2)+((Zw(3))^2))^(3/2))];
    Z1(:,n+1)=ZZ+((h/6)*(k1+(2*k2)+(2*k3)+k4));
    x1(n+1)=Z1(1,n+1);
    x2(n+1)=Z1(2,n+1);
    y1(n+1)=Z1(3,n+1);
    y2(n+1)=Z1(4,n+1);
end

figure(5)
clf
plot(x1,y1,'b')


4 Resultado e interpretación de las gráficas obtenidas