[Home] [Lexikon] [Aufgaben] [Tests] [Kurse] [Begleitmaterial] [Hinweise] [Mitwirkende] [Publikationen]
	Mathematik-Online-Aufgabensammlung: Lösung zu
	Interaktive Aufgabe 1182: N-Körper-Problem

function nbody(p,v,c,T)
% n-body problem 

% Defaults
if ~exist('c','var') | isempty(c)
  c = [100 100 1];
end
if ~exist('T','var') | isempty(T)
  T=6;
end
if ~exist('p','var') | isempty(p)
  p=[ 4,-4,1;0,0,0;0,0,0];
end
if ~exist('v','var') | isempty(v)
  v=[ 0,0,0;-1,1,0;0,0,4];
end

% u=[p_1;v_1;p_2;v_2;...]
u=[p;v];
u=u(:);

% solution and plot
[t,pv] = ode45(@(t,u)NBodyFct(t,u,c),[0,T],u,odeset('reltol',1e-7));

n = length(c); ind = 1 + 6*[0:n-1];
plot3(pv(:,ind),pv(:,ind+1),pv(:,ind+2));

norm(pv(end,end-2:end))

function y = NBodyFct(t,x,c)

% n-body problem 
% p_j'' = sum_{k~=j} c_k (p_k-p_j) / d_jk^3
% x = [p_1; p_1'; p_2; p_2'; ...]

n = length(x)/6; 
y = 0*x;
J = [1:3]; 
for j=1:n; 
   y(J) = x(J+3);
   K = [1:3];
   y(J+3) = zeros(3,1);
   for k=1:n;
      if k~=j; 
         pjk = x(K)-x(J);
         djk = norm(pjk);
         y(J+3) = y(J+3) + c(k)*pjk/(djk^3);
      end;
      K = K+6;
   end;
   J = J+6;
end;

>>nbody
ans =
    9.9802