Runge-kutta results do not align with ODE45 solver, what am I doing wrong?

I am working on simulation the movement of three bodies with mass M, position P, velocity V, acceleration A, and force applied in it F.

I am given the initial values of position and velocity to work with, and according to Newton laws I can acquire the force and the right hand side of equation to implement Runge-Kutta and acquire the position of each body, though when I implement runge-kutta and compare the results to those I get from ODE45, I find them disimilar. What am I doing wrong?

I have attached the file of the instruction I was given and have followed.

The only two files that need to be run is ODE.m and testingrunge.m.

[Jan]

Just some hints for a simplification:

% Replace
P1=[X(1),X(2),X(3)];
% by
P1 = X(1:3);

% Replace
Y_matrix=cell2mat([{X(4),X(5),X(6)};
        {A1(1),A1(2),A1(3)};
        {X(10),X(11),X(12)};
        {A2(1),A2(2),A2(3)};
        {X(16),X(17),X(18)};
        {A3(1),A3(2),A3(3)}]);
Y=(Y_matrix(:));
% by
Y = [X(4:6), A1, X(10:12), A2, X(16:18), A3].';

% Replace
F=((-1)*(M*M1).*(P-P1)/((norm(P-P1).^3)))-(1*M*M2.*(P-P2)./((norm(P-P2)).^3));
% by:
F= -M*M1 * (P-P1) / norm(P-P1).^3 - M*M2 * (P-P2) / (norm(P-P2).^3);

Answer 1

[Jan]

You are using wrong initial values in the Runge Kutta method:

% X(:, 1) = RHS(0, V0);  Nope!
X(:, 1) = V0;

The inital values are the initial values, not the right hand side of them.

I simplified the code to make it easier to debug:

function main
p  = 0.05;
P1 = [0;0;0]; V1 = [p;p;p];
P2 = [1;0;0]; V2 = [-p;p;p];
P3 = [0;0;1]; V3 = [-p;-p;-p];
V0 = [P1; V1; P2; V2; P3; V3];
h  = 0.01;
t0 = 0;
tEnd = 1;

[t, X] = ode45(@RHS, t0:h:tEnd, V0);
figure;
plot(t, X);

[t2, X2] = rungekutta(h, t0, tEnd, V0);
figure;
plot(t2, X2);

end

function Y = RHS(t, X)
    M1 = 1; M2 = 2; M3 = 0.5; 
    
    P1 = X(1:3);
    P2 = X(7:9);
    P3 = X(13:15);
    A1 = Force(P1, M1, P2, M2, P3, M3) / M1;
    A2 = Force(P2, M2, P1, M1, P3, M3) / M2;
    A3 = Force(P3, M3, P1, M1, P2, M2) / M3;
    Y  = [X(4:6); A1; X(10:12); A2; X(16:18); A3];
end

function F = Force(P, M, P1, M1, P2, M2)
    F = -M * M1 * (P-P1) / norm(P-P1).^3 ...
        -M * M2 * (P-P2) / norm(P-P2).^3;
end

function [t, X] = rungekutta(h, t0, Tend, V0)
t  = t0:h:Tend;
X  = zeros(numel(V0), length(t));  % Pre-allocation
X(:, 1) = V0;  % Fixed bug here!
for i = 1:length(t)-1
   k_1 = h * RHS(t(i),         X(:, i));
   k_2 = h * RHS(t(i) + 0.5*h, X(:, i) + 0.5 * h * k_1);
   k_3 = h * RHS(t(i) + h,     X(:, i) - k_1 + 2 * k_2);
   X(:, i+1) = X(:, i) +  (k_1 + 4 * k_2 + k_3) / 6;
end
end

[Nour Butrus]

Thank you so much! I can't believe I was losing my mind for one day over this mistake.

Thank you a lot for simplifying the code!