Runge-kutta results do not align with ODE45 solver, what am I doing wrong?
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaire
Jan
le 17 Fév 2022
Modifié(e) : Jan
le 17 Fév 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);
Réponse acceptée
Jan
le 17 Fév 2022
Modifié(e) : Jan
le 17 Fév 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
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!