help me pls about N-body CODE

13 Ansichten (letzte 30 Tage)
RATCHATA UTAMA
RATCHATA UTAMA am 27 Feb. 2018
Kommentiert: Firas Sheikh am 1 Sep. 2020
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 Kommentare
Deven Mhadgut
Deven Mhadgut am 4 Dez. 2018
Can you share the correct code? I am getting similar errors related to the use of ode45
Firas Sheikh
Firas Sheikh am 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!

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Walter Roberson
Walter Roberson am 27 Feb. 2018

Weitere Antworten (1)

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

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Gravitation, Cosmology & Astrophysics finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by