Using a Maple, simulated a system with three masses and their movement, caused by gravitational forces in 3D space. Two bodies are attracted by a force with a magnitude G x m1 x m2 / r^2 (Newton's Law of Universal Gravitation), where G is the gravitational constant, m1 and m2 are masses of two bodies, and r is the distance between the two. The direction of the force is along the line segment joining two bodies. By equating the equation with Newton's second law (F = ma), we can calculate the acceleration (a) of the body, which is equal to second order derivative of the position of the mass. By solving the deferential equation, we can calculate the position of the bodies at a given time, and plot them in 3 dimensional space. The above figure is generated by running the following maple commands (in red). We can also animate the result by running additional commands in blue at the end.
> d3 := proc(i,j) ((x[i](t)-x[j](t))^2 + (y[i](t)-y[j](t))^2 + (z[i](t)-z[j](t))^2)^(1/2) end:
> v3 := [x,y,z]:
> v3[2][1](t):
> N := 3:
> iN := [seq(i,i=1..N)]:
>
> cde := 0:
for i to N do
i3m := subsop(i=NULL,iN);
for k to 3 do
cde := cde+1;
de[cde] := m[i]*diff(v3[k][i](t),t$2) = add(G*m[i]*m[j]/d3(i,j)^2 * (v3[k][j](t)-v3[k][i](t))/d3(i,j),j=i3m);
od;
od:
> for i to cde do de[i]; od:
>
> des := convert(de,set):
> ic[1] := x[1](0)=0:
> ic[2] := y[1](0)=0:
> ic[3] := z[1](0)=0.5:
> ic[4] := x[2](0)=5:
> ic[5] := y[2](0)=0:
> ic[6] := z[2](0)=0:
> ic[7] := x[3](0)=-2:
> ic[8] := y[3](0)=0:
> ic[9] := z[3](0)=0:
> ic[10] := D(x[1])(0)=.5:
> ic[11] := D(y[1])(0)=0:
> ic[12] := D(z[1])(0)=0:
> ic[13] := D(x[2])(0)=0:
> ic[14] := D(y[2])(0)=.4:
> ic[15] := D(z[2])(0)=0:
> ic[16] := D(x[3])(0)=0:
> ic[17] := D(y[3])(0)=-0.16:
> ic[18] := D(z[3])(0)=0:
> pars := [m[1]=.1,m[2]=3,m[3]=5,G=1]:
>
> ics := convert(ic,set):
> va := [seq(seq(v3[i][j](t),i=1..N),j=1..3)]:
> ds := dsolve(subs(pars,des) union ics,va,numeric):
> pt[0] := ds(0):
> n := 1600:
> for i to n do pt[i] := ds(i*.05); od:
> orb1 := [seq(subs(pt[i],[x[1](t),y[1](t),z[1](t)]),i=0..n)]:
> p1 := plots[pointplot3d](orb1,connect=true,symbol=point,colour=red):
> orb2 := [seq(subs(pt[i],[x[2](t),y[2](t),z[2](t)]),i=0..n)]:
> p2 := plots[pointplot3d](orb2,connect=true,symbol=point,colour=blue):
> orb3 := [seq(subs(pt[i],[x[3](t),y[3](t),z[3](t)]),i=0..n)]:
> p3 := plots[pointplot3d](orb3,connect=true,symbol=point,colour=green):
> plots[display3d]({p1,p2,p3},scaling=constrained,axes=none);
>
> re := COLOUR(RGB,1.,0.,0.):
> bl := COLOUR(RGB,0.,0.,1.):
> gr := COLOUR(RGB,0.,1.,0.):
>
> cir := SYMBOL(CIRCLE):
> for i from 0 to n do
pti := pt[i];
> fr[i] := [POINTS(subs(pti,[x[1](t),y[1](t),z[1](t)]),cir,re),
POINTS(subs(pti,[x[2](t),y[2](t),z[2](t)]),cir,bl),
POINTS(subs(pti,[x[3](t),y[3](t),z[3](t)]),cir,gr),
TITLE(cat("t=",sprintf("%05.2f",subs(pti,t))))];
> od:
> n:fr[n]:
>
> sc := SCALING(CONSTRAINED):
> an := ANIMATE(seq(fr[i],i=0..n)):
> ax := AXESSTYLE(NONE):
> pa := PLOT3D(an,sc,ax):
> pa;
