help me pls about N-body CODE

13 visualizaciones (últimos 30 días)
RATCHATA UTAMA
RATCHATA UTAMA el 27 de Feb. de 2018
Comentada: Firas Sheikh el 1 de Sept. de 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 comentarios
Deven Mhadgut
Deven Mhadgut el 4 de Dic. de 2018
Can you share the correct code? I am getting similar errors related to the use of ode45
Firas Sheikh
Firas Sheikh el 1 de Sept. de 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!

Iniciar sesión para comentar.

Respuesta aceptada

Walter Roberson
Walter Roberson el 27 de Feb. de 2018

Más respuestas (1)

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

Iniciar sesión para comentar.

Categorías

Más información sobre Gravitation, Cosmology & Astrophysics en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by