# Issue with N body problem using ode45

9 views (last 30 days)
Alix DUMAS on 7 Sep 2022
Answered: Alix DUMAS on 7 Sep 2022
Hello,
I am working on a N body problem with N = 4 satellites of mass 10^10kg revolving around the Earth over 3,000s.
Their initial positions are X1 = [100; 0; 0];X2 = [0; 100; 0];X3 = [-100; 0; 0]; X4 = [0; -100; 0]; in meters
And their initial velocities are V1 = [0; 0.1; 0]; V2 = [-0.1; 0; 0]; V3 = [0; -0.1; 0]; V4 = [0.1; 0; 0]. in meters per second
I am using ode45 to solve the problem. To do that I use an X vector defined as follow :
X = (x1 x1dot y1 y1dot z1 z1dot ... xN xNdot yN yNdot zN zNdot) where x,y and z are the coordinates for each satellite.
I put my code below.
close all
clear all
clc
% Number of satellites
N = 4;
% Initial conditions
% X format is [x dx y dy] in m and m/s
X0 = zeros(N*6, 1);
X0(1:6,1) = [100 0 0 0.1 0 0];
X0(7:12,1) = [0 -0.1 100 0 0 0];
X0(13:18,1) = [-100 0 0 -0.1 0 0];
X0(19:24,1) = [0 0.1 -100 0 0 0];
% Propagation over 3,000s
tspan = linspace(0,3000,3000); % [s]
% Change of tolerance for integrator
tol = 1e-4;
options = odeset('RelTol',tol);
% Resolution of the equation with ode45 using RK4 method
[t, x] = ode45(@(t, X) ODE_multiple(t, X, N), tspan, X0, options);
And the ODE_multiple function that computes the derivative of the X vector is coded here :
% Function that represents the equation
% X vector = (x xdot y ydot z zdot)* # of
function dXdt = ODE_multiple(t, X, N)
G = 6.674*10^(-11); % [m^3/(kg.s²)]
mE = 5.9722*10^24; % [kg]
mass = 10^10; % [kg]
dXdt = zeros(length(X),1);
for i = 1:N
r = [X(1 + 6*(i-1)) X(3 + 6*(i-1)) X(5 + 6*(i-1))];
n = norm(r);
ddr = mE/n^3;
for j=1:N
if j~=i
rj = [X(1 + 6*(j-1)) X(3 + 6*(j-1)) X(5 + 6*(j-1))];
ri = [X(1 + 6*(i-1)) X(3 + 6*(i-1)) X(5 + 6*(i-1))];
rij = rj - ri;
nij = norm(rij);
ddr = ddr + mass/nij^3;
end
end
ddr = ddr*G;
% the xdot, ydot and zdot for each satellite
dXdt(1 + 6*(i-1)) = X(2 + 6*(i-1));
dXdt(3 + 6*(i-1)) = X(4 + 6*(i-1));
dXdt(5 + 6*(i-1)) = X(6 + 6*(i-1));
% the x, y and z for each satellite
dXdt(2 + 6*(i-1)) = ddr*X(1 + 6*(i-1));
dXdt(4 + 6*(i-1)) = ddr*X(3 + 6*(i-1));
dXdt(6 + 6*(i-1)) = ddr*X(5 + 6*(i-1));
end
end
The x given by ode45 is a matrix with each positions and velocities as colomns evolving through time (rows). And from what I see in my x variable (see image below the dark grey columns), my velocities are not evolving which is wrong and I don't understand why.
Could you help me with that issue ?

Alix DUMAS on 7 Sep 2022
I ended up finding what was wrong in my code thanks to them :
• First the Earth isn't taken into account here so maybe the bodies are BIG so m = 10^10kg is fine but I had to take off the interaction with the Earth in my ODE function
• Second when I computed my equation to calculate the accelarations I multiplied everything by ri instead of multiplying each term representing body j by the vector rij, so it did not work
So in the end the problem was just to see the interaction of N bodies that started their motion close to each other and then interracted with each other. Thank you for your help !

Steven Lord on 7 Sep 2022
Are you certain on those satellite masses? Looking at Wikipedia the mass of the Great Pyramid of Giza is 6e9 kg. Are these satellites small asteroids about twice the size of the Great Pyramid? Or maybe you're using the wrong units? 1e10 grams would be 1e7 kilograms, or about 5 space shuttles or 1/5 of a Titanic.
They're likely to have a lot of inertia and so accelerate very s l o w l y.
To check this try displaying your data using a longer variable format. If you have this variable open in the Variable Editor, click on the View tab on the Toolstrip and change the number display format. My guess is that these large chunks of rock and/or metal are moving, but not at all quickly.

James Tursa on 7 Sep 2022
Edited: James Tursa on 7 Sep 2022
You need to rethink your problem statement and initial conditions. The initial positions can't possibly be in the neighborhood of 100 meters since that is practically at the center of the Earth. And initial velocities in the neighborhood of 0.1 m/s don't make sense either because that would effectively produce nearly rectilinear motion with everything falling almost straight down. Are these supposed to be "relative" positions and velocities about some reference circular orbit at some altitude?
Also not sure why you are running an N-Body problem in the first place instead of N separate 2-Body problems. Are you trying to simulate the effects of the weak gravitational interaction between 4 small bodies as they orbit in a close group above a much larger Earth? Or ...?
You list the mass of each small body as 10^10 kg, but if I did my calculations correctly a spherical rock with that mass would have a radius of about 98 meters. That means these small bodies are practically touching each other. Is that the intent?
Your current description with the initial positions, velocities, and masses simply does't make sense and you need to correct this and post an update to your question.
SIDE NOTE: It is hard to follow your indexing, but at first glance the velocity derivatives don't look correct. I would have expected to see cumulative vector additions to calculate these derivatives, but I don't see that. I might suggest that you reshape your state vector upon entry to be 3x2N matrix to make your indexing easier. E.g.,
Y = reshape(X,3,[]);
Then the first body position is simply Y(:,1) and its velocity is Y(:,2). The second body Pos & Vel would be Y(:,3) and Y(:,4). The n'th body Pos & Vel is Y(:,2*n-1) and Y(:,2*n). Etc. This will be much easier to write code downstream than what you currently have.
Alternatively you could reshape to a 6xN matrix, e.g.,
Y = reshape(X,6,[]);
Then the first body position is Y(1:3,1) and the velocity is Y(4:6,1). The second body Pos & Vel is Y(1:3,2) and Y(4:6,2). The n'th body Pos & Vel is Y(1:3,n) and Y(4:6,n). Etc.
I suppose you could even get cute and create two helper indexing variables:
pos = 1:3;
vel = 4:6;
Then the n'th body position and velocity would simply be Y(pos,n) and Y(vel,n) which would make your downstream code even more readable.
E.g., using this last scheme you can get all the position derivatives in one line:
dYdt(pos,:) = Y(vel,:);
Then run your loops to fill in the dYdt(vel,:) parts.
After calculating the 3x2N or 6xN matrix derivative dYdt, then reshape it to a column vector for return purposes:
dXdt = dYdt(:);

### Categories

Find more on Ordinary Differential Equations 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!

Translated by