help me pls about N-body CODE

13 views (last 30 days)
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
  2 Comments
Deven Mhadgut
Deven Mhadgut on 4 Dec 2018
Can you share the correct code? I am getting similar errors related to the use of ode45
Firas Sheikh
Firas Sheikh on 1 Sep 2020
I, too, would also like to know if you could share the correct code. I am having trouble with creating my own numerical integration code with ode45. Any help would be appreciated!

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 27 Feb 2018

More Answers (1)

Marco Santambrogio
Marco Santambrogio on 12 Mar 2018
[t,x] = ode45(@(t,x) nBody(t, x, m, G), time, x0);
this is the solution!
  2 Comments
Deven Mhadgut
Deven Mhadgut on 4 Dec 2018
I tried this but its still not working.
Walter Roberson
Walter Roberson on 4 Dec 2018
Deven, what error are you getting, exactly? Please copy it all, everything in red.

Sign in to comment.

Categories

Find more on Gravitation, Cosmology & Astrophysics in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!