MATLAB: Help me pls about N-body CODE

nbodyphysics

nBody.m
function [dx] = nBody(t,x,m,G)
% This is the function that the ode solver needs to solve the n-body problem
% Setting up things
n = length(m); % The number of bodies
dx = zeros(6*n,1); % Holder for the dx
half = 3*n;
% Creating the r matrix
% The idea employed here is that precomputing the repeated values of the
% radius magnitudes cubed saves computation later - and additionally the
% symmetry of the solution may be taken of advantage to reduce the number
% of computations used as well.
r = zeros(n,n);
for i = 1:n
for j = 1:n
if i == j (or) r(i,j) ~= 0 %symbol or not show in this ask.
continue; % The radius between a body and itself is 0
else
r(i,j) = (((x((3*(j-1))+1)-x((3*(i-1))+1)).^2+(x((3*(j-1))+2)-x((3*(i-1))+2)).^2+(x(3*j)-x(3*i)).^2).^(3/2));
r(j,i) = r(i,j);
end
end
end
% Filling in the velocities
for i = 1:half
dx(i) = x(half+i);
end
% Building the accelerations
for i = 1:n
temp = 0; % resetting the holder
% building the x-component
for j = 1:n
if j == i
continue; % The body has no gravitational effect on itself

else
temp = temp + (m(j)/(r(i,j)))*(x((3*(j-1))+1)-x((3*(i-1))+1));
end
end
temp = G*temp;% adjusting for the gravitational constant

dx(half+(3*(i-1))+1) = temp;% setting the dx


temp = 0; % resetting the temp holder

% building the y-component
for j = 1:n
if j == i
continue; %
%The body has no gravitational effect on itself
else
temp = temp + (m(j)/(r(i,j)))*(x((3*(j-1))+2)-x((3*(i-1))+2));
end
end
temp = G*temp; % adjusting for the gravitational constant
dx(half+(3*(i-1))+2) = temp; % setting the dx
temp = 0; % resetting the temp holder
%# building the z-component
for j = 1:n
if j == i
continue; % The body has no gravitational effect on itself
else
temp = temp + (m(j)/(r(i,j)))*(x(3*j)-x(3*i));
end
end
temp = G*temp; %# adjusting for the gravitational constant
dx(half+(3*i)) = temp; % setting the dx
end
end
Twobody.m
% This is a script that simulates two bodies of the same mass orbiting each
% other.
format long;
G = 6.673e-11; % Universal gravitation constant
m = [2, 2]; % Setting both masses as 2 kg
time = [0 1087763]; % running for 1 period worth of seconds
x0 = [-1;0;0;1;0;0;0;-5.775e-6;0;0;5.775e-6;0];
% Starting positions and velocities
[t,x] = ode45(@nBody,time,x0,m,G);
figure();
hold on;
plot(x(:,1),x(:,2));
plot(x(:,4),x(:,5),'*');
Error CODE :
Error using nBody (line 39) Not enough input arguments.
Error in odearguments (line 87) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in TWOBODY (line 9) [t,x] = ode45(@nBody,time,x0,m,G);
I can't run this code.
I'll applied this code for my project.
thank you answer for me 🙂
I copied from Numerical Integration of the n-Body Problem Harrison Brown∗ University of Colorado Boulder, Boulder, Colorado, 80309, USA
Related Question