Matlab - Barridos con matrices de por medio

 
Vista:

Barridos con matrices de por medio

Publicado por Dez (2 intervenciones) el 04/04/2013 22:19:36
Hola buenas.
Estoy intentando hacer un programa para predecir la órbita de asteroides con Euler. Consiste en darle una aceleración, velocidad y posiciones iniciales en x e y, usando el siguiente bucle:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
t(1)=tini;
pos(:,1)=[x0;y0];
a(:,1)=[a0x;a0y];
v(:,1)=[vx0;vy0];
n=1;
 
while ~(pos(1,n)>0 && pos(2,n)<0 )
    n=n+1;
    t(n)=t(n-1)+da;
 
    a(:,n)=-G*M.*pos(:,n-1)/(sum(pos(:,n-1).^2)).^(3/2);
    v(:,n)=v(:,n-1)+a(:,n)*da;
    pos(:,n)=pos(:,n-1)+(v(:,n-1)+v(:,n))/2*da;
 
end


Los valores se guardan en matrices 2xn (en la primera fila los valores de x y en la segunda los de y)
El problema llega cuando quiero meter la velocidad como un vector e intentar hacer un bucle que haga lo mismo que el anterior pero con un barrido de todas las velocidades. He intentado meter un "for" que las recorra, pero no puedo poner v(i,:,n-1) por ejemplo. ¿Cómo podría solucionarlo? Gracias
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder
Imágen de perfil de JOSE JEREMIAS CABALLERO
Val: 6.975
Oro
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

Barridos con matrices de por medio

Publicado por JOSE JEREMIAS CABALLERO (5917 intervenciones) el 06/04/2013 01:05:57
Me hubiera gustado ejecutar tu código, en tu código que pones , no mencionas los valores numéricos de varios de tu datos y eso ocasiona no pueda ejecutar y menos entender tu código. Y va ser difícil que alguien te ayuda con la informacion q estas dando, los errores pueden ser como por ejemplo, distinto tamaño de vector pos y el vector v. (creo que son matrices, hay muchos vacíos en tu pregunta que haces).
Conclusión: mejorar la pregunta.

Saludos.
JOSE JEREMÍAS CABALLERO
Asesorías en Matlab
programador en matlab
Servicios de programación matlab
[email protected]

http://matlabcaballero.blogspot.com

http://www.lawebdelprogramador.com/foros/Matlab/1371532-FORMA_DE_APRENDER_MATLAB.html
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Barridos con matrices de por medio

Publicado por Dez (2 intervenciones) el 06/04/2013 14:03:17
Tienes razón, perdona, aquí va el programa completo.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
home
clear all
close all
 
M=5.98e24;
m=500000;
G=6.67e-11;
rt=3670000;
 
x0=-200*rt;
y0=0;
r0=sqrt(x0^2+y0^2);
vector0=[x0;y0];
 
ango0=0*(pi/180);
da=10;
 
c=(-G*M)/(r0^3);
a0x=c.*(x0);
a0y=c.*(y0);
a0=c.*(vector0);
 
vx0=5000*(10/36);
vy0=1*(10/36);
v0=(sqrt(vx0^2+vy0^2));
 
tini=0;
tfin=200;
npasos=(tfin-tini)/da+1;
 
%%
 
 
 
t(1)=tini;
pos(:,1)=[x0;y0];
a(:,1)=[a0x;a0y];
v(:,1)=[vx0;vy0];
n=1;
 
while ~(pos(1,n)>0 && pos(2,n)<0 )
    n=n+1;
    t(n)=t(n-1)+da;
 
    a(:,n)=-G*M.*pos(:,n-1)/(sum(pos(:,n-1).^2)).^(3/2);
    v(:,n)=v(:,n-1)+a(:,n)*da;
    pos(:,n)=pos(:,n-1)+(v(:,n-1)+v(:,n))/2*da;
 
 
end
 
choque=find( (sum(pos(:,n-1).^2)).^(1/2) <= rt );
posc=pos(:,choque);
 
 
plot(pos(1,:),pos(2,:),'k','LineWidth',2)
hold on
plot(0,0,'r.','MarkerSize',20)
plot(x0,y0,'b.','MarkerSize',20)
 
axis equal


En este caso es para una única velocidad. Mi problema llega al convertir las velocidades en vectores (con distintas velocidades)
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar
Imágen de perfil de JOSE JEREMIAS CABALLERO
Val: 6.975
Oro
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

Barridos con matrices de por medio

Publicado por JOSE JEREMIAS CABALLERO (5917 intervenciones) el 06/04/2013 16:33:14
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
home
clear all
close all
 
M=5.98e24;
m=500000;
G=6.67e-11;
rt=3670000;
 
x0=-200*rt;
y0=0;
r0=sqrt(x0^2+y0^2);
vector0=[x0;y0];
 
ango0=0*(pi/180);
da=10;
 
c=(-G*M)/(r0^3);
a0x=c.*(x0);
a0y=c.*(y0);
a0=c.*(vector0);
 
vx0=5000*(10/36);
vy0=1*(10/36);
v0=(sqrt(vx0^2+vy0^2));
 
tini=0;
tfin=200;
npasos=(tfin-tini)/da+1;
 
%%
 
 
 
t(1)=tini;
pos(:,1)=[x0;y0];
a(:,1)=[a0x;a0y];
v(:,1)=[vx0;vy0];
n=1;
 
while ~(pos(1,n)>0 && pos(2,n)<0 )
    n=n+1;
    t(n)=t(n-1)+da;
 
    a(:,n)=-G*M.*pos(:,n-1)/(sum(pos(:,n-1).^2)).^(3/2);
    v(:,n)=v(:,n-1)+a(:,n)*da;
    pos(:,n)=pos(:,n-1)+(v(:,n-1)+v(:,n))/2*da;
 
 
end
 
choque=find( (sum(pos(:,n-1).^2)).^(1/2) <= rt );
posc=pos(:,choque);
 
 
plot(pos(1,:),pos(2,:),'k','LineWidth',2)
hold on
plot(0,0,'r.','MarkerSize',20)
plot(x0,y0,'b.','MarkerSize',20)
 
axis equal



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
EJECUCION
>> asteroide_de_euler
  Name         Size                  Bytes  Class     Attributes
 
  G            1x1                       8  double
  M            1x1                       8  double
  a            2x39083              625328  double
  a0           2x1                      16  double
  a0x          1x1                       8  double
  a0y          1x1                       8  double
  ango0        1x1                       8  double
  c            1x1                       8  double
  choque       1x1                       8  double
  da           1x1                       8  double
  m            1x1                       8  double
  n            1x1                       8  double
  npasos       1x1                       8  double
  pos          2x39083              625328  double
  posc         2x1                      16  double
  r0           1x1                       8  double
  rt           1x1                       8  double
  t            1x39083              312664  double
  tfin         1x1                       8  double
  tini         1x1                       8  double
  v            2x1x39083            625328  double
  v0           1x1                       8  double
  vector0      2x1                      16  double
  vx0          1x1                       8  double
  vy0          1x1                       8  double
  x0           1x1                       8  double
  y0           1x1                       8  double



Saludos.
JOSE JEREMÍAS CABALLERO
Asesorías en Matlab
programador en matlab
Servicios de programación matlab
[email protected]

http://matlabcaballero.blogspot.com

http://www.lawebdelprogramador.com/foros/Matlab/1371532-FORMA_DE_APRENDER_MATLAB.html
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar