Problema restringido de los tres cuerpos

Se usa el método de Runge-Kutta de cuarto orden (RK4), para solucionar el problema junto con Máxima. Se asume que la orbita de los dos cuerpos masivos es una circunferencia y se usa la proporciones entre las masas de la Tierra y la Luna en el ejemplo, refiriéndonos como 1 a la masa de la Tierra y a la distancia de la Tierra a la Luna. Debido a estos ajustes las unidades que se manejan en las ecuaciones cambian y hay que tenerlo en cuenta.

El problema se resuelve usando Maxima, el video es simplemente una unión de las imágenes sucesivas producidas por el programa. Las instrucciones para ese propósito no están en el código mostrado.

Como un uso practico, el código puede ser usado para probar las orbitas estables e inestables para el problema.

/* [wxMaxima: input   start ] */

/*
    Masa de la Tierra igual a 1.
    Tiempo medido en periodos.
    Distancia Tierra-Luna igual a 1.
    M = masa de la Luna.
    dt = espacio entre intervalo temporal.
    iter = numero de iteraciones base k.
    x0, y0, v0x, v0y = pocicion y velocidad inicial.
*/

M:1/81.3$n:sqrt(M+1)$
rug4(dt,iter,x0,y0,v0x,v0y):=(
    kill(grap),kill(grap2),

    iter:iter*1000,

    axf(tx,ty):=bfloat(-tx/(tx^2+ty^2)^(3/2)+M*((xl-tx)/(sqrt((tx-xl)^2+(ty-yl)^2))^3-xl)),
    ayf(tx,ty):=bfloat(-ty/(tx^2+ty^2)^(3/2)+M*((yl-ty)/(sqrt((tx-xl)^2+(ty-yl)^2))^3-yl)),

    for i:0 thru iter do(

        xl:cos(n*i*dt),
        yl:sin(n*i*dt),

        if mod(i,33)=0 then (grap[i+1]:[x0,y0],grap2[i+1]:[xl,yl]),

        kx1:dt*v0x,
        ky1:dt*v0y,

        lx1:dt*axf(x0,y0),
        ly1:dt*ayf(x0,y0),

        kx2:dt*(v0x+lx1/2),
        ky2:dt*(v0y+ly1/2),

        lx2:dt*axf(x0+kx1/2,y0+ky1/2),
        ly2:dt*ayf(x0+kx1/2,y0+ky1/2),

        kx3:dt*(v0x+lx2/2),
        ky3:dt*(v0y+ly2/2),

        lx3:dt*axf(x0+kx2/2,y0+ky2/2),
        ly3:dt*ayf(x0+kx2/2,y0+ky2/2),

        kx4:dt*(v0x+lx3),
        ky4:dt*(v0y+ly3),

        lx4:dt*axf(x0+kx3,y0+ky3),
        ly4:dt*ayf(x0+kx3,y0+ky3),

        x0:bfloat(x0+(kx1+2*kx2+2*kx3+kx4)/6),
        v0x:bfloat(v0x+(lx1+2*lx2+2*lx3+lx4)/6),

        y0:bfloat(y0+(ky1+2*ky2+2*ky3+ky4)/6),
        v0y:bfloat(v0y+(ly1+2*ly2+2*ly3+ly4)/6)
    ),
    print("N° de periodos: ",float(dt*iter*n/(2*%pi))),

    plot2d(
        [[discrete,listarray(grap)],[discrete,listarray(grap2)]],
        [x,-1.2,1.2],[y,-1.2,1.2],
        [legend,"Objeto de prueba","Luna"]
    )

)$

/* [wxMaxima: input   end   ] */

Video en YouTube