Runge-kutta results do not align with ODE45 solver, what am I doing wrong?
3 views (last 30 days)
Show older comments
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 Comment
Jan
on 17 Feb 2022
Edited: Jan
on 17 Feb 2022
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);
Accepted Answer
Jan
on 17 Feb 2022
Edited: Jan
on 17 Feb 2022
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
More Answers (0)
See Also
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!