Image:N-body problem (3).gif

[edit] Summary

The chaotic movement of 3 interacting particles.

This is an animation of the (not-reduced) three-body problem. The center of view is the mass center of the three particles.

Calculated numerically with Maple 10. Note: the 'n' parameter (currently '3') can be adapted at will.

> restart;with(RandomTools);with(plottools):with(plots);
> n:=8;
> e:={}:for i from 1 to n do
   s:=[]:for c from 1 to n do;  if (not(c=i)) then
   s:=[op(s),G*mass[i]*mass[c]*(x[c](t)-x[i](t))/sqrt(((x[c](t)- x[i](t))^2+(y[c](t)-y[i](t))^2+(z[c](t)-z[i](t))^2))^3]: end if od:
   e:={op(e),mass[i]*diff(x[i](t),t$2)=add(s[v],v=1..n-1)}:
  
   s:=[]:for c from 1 to n do;  if (not(c=i)) then
   s:=[op(s),G*mass[i]*mass[c]*(y[c](t)-y[i](t))/sqrt(((x[c](t)- x[i](t))^2+(y[c](t)-y[i](t))^2+(z[c](t)-z[i](t))^2))^3]: end if od:
   e:={op(e),mass[i]*diff(y[i](t),t$2)=add(s[v],v=1..n-1)}:
  
   s:=[]:for c from 1 to n do;  if (not(c=i)) then
   s:=[op(s),G*mass[i]*mass[c]*(z[c](t)-z[i](t))/sqrt(((x[c](t)- x[i](t))^2+(y[c](t)-y[i](t))^2+(z[c](t)-z[i](t))^2))^3]: end if od:
   e:={op(e),mass[i]*diff(z[i](t),t$2)=add(s[v],v=1..n-1)}:
  od:e;
> G:=200;
> SetState(state=12);
> pos:=[]:for i from 1 to n do pos:=[op(pos),Generate(list(integer(range=-10..10), 3))] od;
> vel:=[]:for i from 1 to n do vel:=[op(vel),Generate(list(integer(range=-5..5), 3))] od;
> mass:=Generate(list(integer(range=5..35), n));
> with(LinearAlgebra):
> beginvwn:=NULL:
   for i from 1 to n do beginvwn:=beginvwn,x[i](0)=pos[i][1],y[i](0)=pos[i][2],z[i](0)=pos[i][3] od:
   for i from 1 to n do beginvwn:=beginvwn,D(x[i])(0)=vel[i][1],D(y[i])(0)=vel[i][2],D(z[i])(0)=vel[i][3] od:
  beginvwn;
> var:=NULL:for i from 1 to n do var:=var,x[i](t),y[i](t),z[i](t) od:var;
> astappen:=500:timescale:=0.02:
  opl:=dsolve([op(e),beginvwn],numeric,output=array([seq(m*timescale,m=0..astappen-1)]),maxfun=500000):
> sl:=convert(opl[1,1],list);
> dp:=ListTools[Flatten]([seq([2+2*k,2+2*k+2*n,2+2*k+4*n],k=0..n-1)]);
> data:=convert(opl[2,1],listlist):
> mc:=[seq(map(list->sum(list[dp[cd+3*cs]]*mass[cs+1],cs=0..n-1)/sum(mass[cs],cs=1..n),data),cd=1..3)]:
> aschijfjes:=astappen;p:=[]:
  for k from 1 to aschijfjes do
  pts:=NULL:for i from 1 to n do pts:=pts,point([seq(data[k][dp[c+3*(i-1)]],c=1..3)],symbolsize=mass[i],color=red,symbol=DIAMOND) od:
  pp:=display(pts);
  p:=[op(p),pp];
  end do:
> display(p,insequence=true,scaling=constrained,axes=boxed);
> aschijfjes:=astappen;p:=[]:
  for k from 1 to aschijfjes do
  pts:=NULL:for i from 1 to n do pts:=pts,point([seq(mc[c][k]-data[k][dp[c+3*(i-1)]],c=1..3)],symbolsize=mass[i],color=red,symbol=DIAMOND) od:
  pp:=display(pts);
  p:=[op(p),pp];
  end do:
> display(p,insequence=true,scaling=constrained,axes=boxed);

[edit] Licensing

I, the copyright holder of this work, hereby release it into the public domain. This applies worldwide. In case this is not legally possible: I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law. Afrikaans | Alemannisch | Aragonés | العربية | Asturianu | Български | Català | Česky | Cymraeg | Dansk | Deutsch | Eʋegbe | Ελληνικά | English | Español | Esperanto | Euskara | Estremeñu | فارسی | Français | Galego | 한국어 | हिन्दी | Hrvatski | Ido | Bahasa Indonesia | Íslenska | Italiano | עברית | Kurdî / كوردی | Latina | Lietuvių | Latviešu | Magyar | Македонски | Bahasa Melayu | Nederlands | ‪Norsk (bokmål)‬ | ‪Norsk (nynorsk)‬ | 日本語 | Polski | Português | Ripoarisch | Română | Русский | Shqip | Slovenčina | Slovenščina | Српски / Srpski | Svenska | ไทย | Tagalog | Türkçe | Українська | Tiếng Việt | Walon | ‪中文(简体)‬ | ‪中文(繁體)‬ | zh-yue-hant | +/-

• N-body problem