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

Question

Nour Butrus 2022-2-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1652060-runge-kutta-results-do-not-align-with-ode45-solver-what-am-i-doing-wrong

评论： Nour Butrus 2022-2-17

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2022-2-17

编辑：Jan 2022-2-17

在 MATLAB Online 中打开

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 2022-2-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1652060-runge-kutta-results-do-not-align-with-ode45-solver-what-am-i-doing-wrong#answer_898420

编辑：Jan 2022-2-17

在 MATLAB Online 中打开

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