Page 37 - Tlahuizcalli CB-28_Neat

cálculo parcial de las y es (sólo se calcula la t(1) = ta;
fuerza entre la masa i-ésima y la masa j-ésima):
xi(1) = ra(1,1);
function fuerzaijvec = fuerza(xi,yi,xj,yj) xi(2) = va(1,1);
xi(3) = ra(2,1);
xi(4) = va(2,1);
rijvec=[xi-xj,yi-yj];
rijvec=rijvec/(sqrt(rijvec(1)^2+rijvec(2)^2))^3; xi(5) = ra(3,1);
fuerzaijvec=rijvec;
xi(6) = va(3,1);
xi(7) = ra(4,1);
endfunction
xi(8) = va(4,1);
xi(9) = ra(5,1);
Nótese que estamos haciendo uso de la xi(10) = va(5,1);
capacidad de operaciones sobre arreglos de xi(11) = ra(6,1);
SCILAB (la rutina regresa en una sola ejecución las xi(12) = va(6,1);
componentes y de la fuerza). Por su parte, la
rutina que evalúa es: x1(:,1) = xi(1);
vx1(:,1) = xi(2);
function f = evalua(t,x) y1(:,1) = xi(3);
vy1(:,1) = xi(4);
G=1; x2(:,1) = xi(5);
m1=1; vx2(:,1) = xi(6);
m2=1; y2(:,1) = xi(7);
m3=1; vy2(:,1) = xi(8);
x3(:,1) = xi(9);
fuerza12vec=-G*m1*m2*fuerza(x(1),x(3),x(5),x(7)); vx3(:,1) = xi(10);
fuerza21vec=-fuerza12vec; y3(:,1) = xi(11);
vy3(:,1) = xi(12);
fuerza13vec=-G*m1*m3*fuerza(x(1),x(3),x(9),x(11));
fuerza31vec=-fuerza13vec; i=1;
h=(tb – ta)/(n-1);
while i < n
fuerza23vec=-G*m2*m3*fuerza(x(5),x(7),x(9),x(11));
fuerza32vec=-fuerza23vec;
k1 = evalua(t(i), xi);
f(1)=x(2);
f(2)=1/m1*(fuerza12vec(1)+fuerza13vec(1)); k2 = evalua(t(i)+h/2, xi+h*k1/2);
f(3)=x(4);
f(4)=1/m1*(fuerza12vec(2)+fuerza13vec(2)); k3 = evalua(t(i)+h/2, xi+h*k2/2);
f(5)=x(6);
f(6)=1/m2*(fuerza21vec(1)+fuerza23vec(1)); k4 = evalua(t(i)+h, xi+h*k3);
f(7)=x(8);
f(8)=1/m2*(fuerza21vec(2)+fuerza23vec(2)); i = i+1;
f(9)=x(10);
f(10)=1/m3*(fuerza31vec(1)+fuerza32vec(1)); t(i) = t(i-1) + h;
f(11)=x(12);
f(12)=1/m3*(fuerza31vec(2)+fuerza32vec(2)); xi = xi + 1/6*h*(k1+2*k2+2*k3+k4);

endfunction x1(i) = xi(1);
vx1(i) = xi(2);
Por último, la rutina principal que implementa el y1(i) = xi(3);
método de RK4 para nuestro problema es la vy1(i) = xi(4);
siguiente: x2(i) = xi(5);
vx2(i) = xi(6);
y2(i) = xi(7);
function [t,x1,vx1,y1,vy1,x2,vx2,y2,vy2,x3,vx3,y3,vy3,n] =
runge2_3(ta,tb,ra,va) vy2(i) = xi(8);
x3(i) = xi(9);
vx3(i) = xi(10);
n = 20001;

36
Año 10 Núm. 28 enero-abril 2024 Tlahuizcalli ISSN: 2448-7260